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ABSTRACT 


Two  companion  algorithms  are  developed  for  constructing  Pade  fractions  along  an  off- 
diagonal  path  of  the  Pade  table  for  a  function  —  A(z)  /  B(z),  where  A(z )  and  B(z)  are  formal 
power  series  over  a  field. 

One  of  the  algorithms  computes  the  first  n  Pade  fractions  along  the  off-diagonal  in  time 
0{n 2).  When  A(z)  and  B(z)  are  finite  power  senes  (i.e.,  polynomials),  it  is  shown  that  the 
algorithm  is  equivalent  to  Euclid’s  extended  algorithm  for  computing  greatest  common  divisors. 

The  other  algorithm,  a  generalization  of  the  first,  proceeds  along  the  off-diagonal  in  qua¬ 
dratic  steps,  and  is  of  complexity  0(n  log2  n).  When  A(z)  and  B(z )  are  polynomials,  the 
second  algorithm  becomes  a  fast  Euclid’s  extended  algorithm  for  computing  greatest  common  divi¬ 
sors.  The  algorithm  is  of  the  same  complexity  as  other  fast  greatest  common  divisor  methods,  but 
its  iterative  nature  provides  a  practical  advantage  during  implementation. 

The  algorithms  may  also  be  used  for  computing  Pade  fractions  along  an  anti-diagonal  path 
of  the  Pade  table,  as  well.  The  fast  algorithm  is  of  the  same  complexity  as  other  fast  algorithms 
for  anti-diagonal  computations.  However,  it  has  the  advantage  of  being  able  to  determine  easily 
any  specific  Pade  fraction  along  the  anti-diagonal. 

Finally,  it  is  shown  that  two  successive  Pade  fractions  can  be  used  to  obtain  the  inverse  of 
Hankel  and  Toeplitz  matrices  in  time  0(n  log2  n). 
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CHAPTER  1 


INTRODUCTION 


The  Pade  table  of  a  fomal  power  series 


A{z)  =  'ZaJ 
i  -  o 


is  a  doubly  infinite  array  of  rational  functions 


Vjj) 

V„  (z) 


2“^ 


i  -  0 

/I 


2v“ 


/-o 


(1.1) 


(1.2) 


determined  in  such  a  manner  that  the  Maclaurin  expansion  of  Umn(z)/Vmn(z)  agrees  with  A{z)  as 
far  as  possible.  The  power  series  A(z)  is  said  to  be  normal  if  ,  for  each  pair  ( m,n ),  tliis  agreement 
is  exact  through  the  power  zm+n.  The  foundation  for  the  development  of  Pade  theory  was  laid  by 
Cauchy (1821)  in  his  famous  "Cour  d’Analyse".  Later,  Frobenius(1881)  developed  the  basic  algo¬ 
rithmic  aspects  of  the  theory,  and  Pade(1892)  treated  in  detail  certain  abnormal  cases. 

Since  Pade’s  time,  Pade  tables  have  become  a  classical  tool  of  analysis.  Their  analytical  pro¬ 
perties  have  been  studied  in  great  depth  and  are  surveyed,  for  example,  by  Gragg  [GRA  74]  and  by 
Baker  [BAK15].  Traditionally,  it  is  assumed  that  the  coefficients  in  (1.1)  and  (1.2)  lie  in  the  field 
of  complex  numbers,  and  that  the  power  series  and  the  rational  functions  are  to  be  evaluated  at 
certain  points  in  the  complex  plane. 

Although  the  results  obtained  in  this  thesis  are  likely  to  have  an  impact  in  an  analytical  (or 
numerical)  setting,  the  effects  of  this  impact  are  not  examined.  The  issues  addressed  are  strictly 
algebraic  ones;  that  is,  no  consideration  is  given  to  the  goodness  of  the  approximation  of  (1.2)  to 
(1.1).  Instead,  the  objective  is  to  provide  an  effective  tool  to  algebraically  manipulate  rational  func- 
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irons  as  truncated  power  series  (for  which  the  cost  of  operations  is  relatively  cheap),  and  to 
transform  back  to  rational  form  on  request.  It  is  assumed  that  the  coefficients  lie  in  an  arbitrary 
field. 

The  coefficients  of  the  rational  function  (1.2)  satisfy  an  under-determined  system  of  linear 
equations,  known  as  the  Hankel  system.  Various  properties  of  the  family  of  solutions  to  the 
Hankel  sytem  are  described  in  Qiapter  2.  All  the  results  given  are  those  of  Pade,  and  are  included 
for  the  sake  of  completeness  and  for  reference  ease.  The  Hankel  system  can  be  solved  directly  for 
the  coefficients  of  Umn (z)  and  V^z).  With  coefficients  over  a  field,  Rissannen  [77S73],  for  exam¬ 
ple,  provides  an  algorithm  which  requires  0(n 2)  arithmetic  operations  over  the  field.  On  the  other 
hand,  with  coefficients  over  an  arbitrary  integral  domain,  Geddes  [G£D79]  gives  a  fraction-free 
algorithm  which  requires  0(n3)  arithmetic  operations  over  the  integral  domain. 

Various  relationships  are  known  to  exist  between  neighboring  elements  in  the  Pade  table. 
These  relationships  have  been  used  to  derive  numerous  0{n2)  methods  for  computing  a  sequence 
of  elements  in  the  Pade  table.  A  survey  and  comparison  of  these  methods  are  given  in  Brezinski 
[BRZ76] ,  Qaessens  [CLA15]  and  Wynn  [WT7V60].  All  these  methods  have  a  major  flaw;  they  may 
fail  in  the  abnormal  case.  In  Qiapter  3,  a  new  relationship  between  elements  lying  along  an  off- 
diagonal  path  in  the  Pade  table  is  derived.  This  leads  to  yet  another  0{n2)  method,  however,  the 
new  method  succeeds  in  the  abnormal  case.  Furthermore,  if  the  coefficient  field  has  an  appropriate 
n-th  root  of  unity  (which  permits  fast  multiplication  and  division  of  polynomials),  the  asymptotic 
complexity  of  the  algorithm  becomes  0(n  log2  n). 

The  new  algorithm  can  be  applied  to  the  quotient  of  two  power  series.  Then,  in  particular, 
it  can  be  applied  to  the  quotient  of  two  finite  power  series  (i.e.  the  quotient  of  two  polynomials). 
In  Qiapter  4,  it  is  shown  that  if  all  elements  along  a  specific  off-diagonal  of  the  Pade  table  are 
computed,  then  the  new  algorithm  is  equivalent  to  Euclid’s  extended  algorithm  for  computing 
greatest  common  divisors.  Furthermore,  if  fast  polynomial  operations  can  be  performed,  the  new 
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algorithm  can  compute  the  greatest  common  divisor  of  two  polynomials  in  0(n  log2  n)  arithmetic 
operations.  The  algorithm  has  three  advantages  over  the  other  fast  methods  (Moenck  [A/0  £73], 
Aho  et  al[A?/074],  and  Brent  et  al  [5/?£80])  for  computing  greatest  common  divisors.  It  is  basi¬ 
cally  an  iterative  algorithm  rather  than  a  recursive  one,  and  consequently,  significant  cost  savings 
can  result  during  implementation.  Secondly,  it  produces  intermediate  polynomial  remainder 
sequences  as  a  by-product,  which  is  a  valuable  feature  for  some  applications.  Finally,  various 
details  about  the  nature  of  its  behavior  are  easier  to  comprehend. 

The  algorithm  can  be  applied  to  the  quotient  of  the  reciprocals  of  two  truncated  power  series 
(polynomials).  It  is  shown  in  Section  4.4  that  this  yields  successive  elements  along  an  anti-diagonal 
path  of  the  Pade  table. 

In  Chapter  5,  it  is  shown  that  two  successive  elements  along  the  diagonal  path  of  the  Pade 
table  can  be  used  to  find  the  inverse  of  a  Hank  el  matrix.  Therefore,  if  fast  polynomial  operations 
are  possible,  the  inverse  of  a  Hankel  matrix  of  order  n  can  be  determined  in  0(n  log2  n)  arith¬ 
metic  operations.  The  new  inversion  method  handles  abnormalities  with  greater  ease,  and  less  cost, 
than  does  Brent  et  al’s  algorithm  [5££80].  Furthermore,  by  proceeding  along  an  anti-diagonal 
path,  the  algorithm  can  be  used  to  compute  the  inverse  of  a  Toeplitz  matrix  of  order  n  in  time 
0(n  log2  n). 


CHAPTER  2 


PADE  THEORY 


2.1.  Introduction 

This  chapter  examines  the  theoretical  background  of  Pade  theory.  None  of  the  results  are 
new;  however,  some  of  the  proofs  may  be  original.  All  of  the  proofs  are  obtained  directly  from  the 
properties  of  Hankel  systems. 

The  highlight  of  this  chapter  is  Corollary  2.6,  due  to  Pade,  which  completely  describes  the 
family  of  solutions  to  the  Hankel  system.  We  define  one  specific  member  of  this  family  to  be  a 
scaled  Pade  fraction.  Scaled  Pade  fractions  exist  uniquely,  and  are  fundamental  to  the  development 
of  subsequent  algorithms. 


2.2.  Pade  Forms 


2.2.1.  Definitions 

The  class  P  of  formal  power  series  over  a  field  F  consists  of  expressions  of  the  form 

A{z)  =  2^ 

i- o 

with  coefficients  a,  e  F  .  We  denote  the  units  of  P  by 

U  =  { A(z)  =  I  ao  *  0  ,  A(z)  e  P  }  . 

<-o 

Associated  with  each  such  unit  is  a  set  of  rational  functions  defined  as  follows: 
Definition:  Let  A(z)  t  U,  and  let  m  and  n  be  non-negative  integers.  The  rational  form 
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_  U0  +  UXZ  +  •  •  •  +  UmZm 
^mn(z)  V0  +  VyZ  +  •  •  ■  +  VnZ" 

(2.1) 

is  called  a  ?ade  form  of  type  (m,  n )  for  A(z)  if 

(a)  Vmn(z)  *  0  ,  and 

(2.2) 

(b)  A(j)-V„(r)  -  Um(z)  =  0(2-"-')  . 

(2.3) 

The  (algebraic)  O-symbol  indicates  that  the  right  side  is  a  power  series  beginning  with  the  power 
zm+n+k  + 1  }  0  ^  jt:  <  oo  ;  k  =  +  *>  means  that  A{z)  •  V^z)  -  U^z)  =  0  . 

From  an  analytical  point  of  view,  this  means  that  the  rational  form  (z)  /  Vmn  (z)  is 

determined  so  that  its  Maclaurin  expansion  agrees  with  A(z)  as  far  as  possible.  That  is,  equation 
(2.3)  corresponds  to  a  linear  system  of  m+n+\  equations  in  the  m+n+2  unknowns 

uo,  v0,  •  •  •  ,  v*  as  follows:1 


n 


2 ai-j‘vj 
j-0 


ui>  i  =  0,  1,  •  ■  •  ,  m, 

0,  i  —  m  +  1,  •  •  •  ,  m  +  n  . 


(2.4) 


Equivalently,  (2.4)  can  be  written  as  the  following  subsystems  for  the  polynomials  Umn{z )  and 

VUz):2 


—  n 

•  •  • 

vn 

• 

= 

■ 

• 

■ 

u0 

a-n 

a0 

v° 

1  By  convention,  at  =  0  if  i  <  0. 

2  Throughout  this  thesis,  the  notation  P(z)  =  p0  +  p^z  +  •  •  •  +  pkzt  for  a  polynomial  P, 
shall  be  used  interchangeably  with  the  vector  notation  P  =  (pk,  pk-U  •  •  ■  ,  p0 )'  for  its  coeffi¬ 
cients. 


6 


On  —  n  + 1 


am 


^ m  + 1 


+n 


V„" 

- 1 

o 

v0 

•  ■  o 

_ 1 

The  matrix 


(2.6) 


Hmn  = 


^ m  -n  +  1 


Om 


“m 


+«  — 1 


(2.7) 


associated  with  the  power  series  A(z)  is  called  a  Hankel  matrix.  With  this  definition,  the  system 
(2.6)  becomes 


'  v„  ' 

^m  +  l 

• 

=  Vg  ■ 

• 

.  V1 

+  n 

(2.8) 


Denote  the  determinant  of  by  det  (//„„,),  and  let  det(//w  0)  =  1. 


2.2.2,  Existence  and  Non-uniqueness 


Theorem  2.1  (Frobenius):  There  always  exists  a  non-trivial  solution  of  the  systems  (2.5)  and 
(2.6).  Equivalently,  Pade  forms  of  type  ( m ,  n )  for  A(z)  e  U  always  exist. 

Proof:  The  existence  of  =  (v„,  •  •  •  ,  v0)  satisfying  (2.6)  follows  immediately  from  the  fact 
that  the  n  +  1  columns  of  the  matrix  in  (2.6)  are  vectors  of  length  n  and  must  therefore  be 
linearly  dependent.  The  vector  U1^  =  («„,,  •  •  •  ,  u0)  satisfying  (2.5)  is  then  obtained  simply  by 
multiplication.  a 


. 


7 


It  is  clear  that  if  U ^  and  V^.  is  a  solution  of  (2.5)  and  (2.6),  then  so  is  a-U ^  and  a  V^, 
where  a  is  a  non-zero  constant.  The  non-uniqueness  of  the  solution,  however,  is  more  profound 
than  this,  and  the  next  few  results  identify  the  true  nature  of  the  family  of  solutions. 


Lemma  2.2:  Let  Umn(z)/Vmn(z)  be  a  Pade  form  of  type  (m,  n )  for  A(z)  e  U,  and  let  D(z)  be  a 
non-zero  polynomial  of  degree  d(D)  <  X,  where  X  =  min{m  —  d(Umn),  n  —  d(V„„,)}.  Then 
U'mn(z)/V'mn(z),  where 

=  Umn(z)D(z)  and  V'^z)  =  Vmfl(z)D(z)  , 
is  also  a  Pade  form  of  type  (m,  n )  for  A(z). 


Proof:  With  the  U^z)  and  Vnu,(z)  defined  above,  equation  (2.5)  becomes 


'  Ox  ■ 

Wjti  —  ft 

'  Ox  • 

Wm-X 

— : 

• 

• 

V*-X 

• 

W0 

a-n 

a0 

Vo 

(2.9) 


where  Ox  refers  to  the  zero  vector  of  length  X .  Multiplication  by  the  polynomial  D  then  results  in 


0 


0 

<*x 

do 


■  Ox  ■ 

Wm-X 

w0 

^771 


a0 


Ox 


V0 
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That  is, 


an  -  n  •  *  *  dm 

dQ  •  •  •  dx 

'Ox  ' 

•  • 

d x 

v„-x 

• 

0 

• 

d-H  '  '  *  a0 

dQ 

v° 

u 


t 

m 


dm  —  * 


Om 


«  0 


a. 


a0 


Therefore,  U'mn  and  V'^  satisfy  equation  (2.5).  The  fact  that  V'mn  also  satisfies  (2.6)  follows 


from 


dm  — n  +  1 


am 


dm  + 1 

v'n 

dm  +n 

v'o 

dm-n  +  l  '  dm+i 

do  ■  "  ■  dy 

°x 

• 

‘  ^ 

1 

>■ 

• 

0 

• 

dm  +  „ 

r - 

v° 

<4  •  •  • 

<4 

dm -«+l 

d m  + 1 

VB-X 

Q 

• 

0 

. 

. 

Vo 

<4 

•  •  • 

do 

dm 

dm  +n 

Ox 

dm -«-X+ 1 

’  dm- 

-  X  + 1 

<4 

Ox  ' 

0 

■ 

• 

• 

• 

dm  -  n 

dn 

v„-x 

• 

• 

dm-n  +  1 

dm  + 1 

• 

0 

<4 

.  .  . 

do 

■ 

• 

V0 

dm 

dm  +n  — 

X 

. 
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0 


0 


(2.10) 


In  the  last  equality  of  (2.10),  use  was  made  both  of  equation  (2.9)  and  of  equation  (2.6)  with 
VL  =  (Ox ,  v„  _  x ,  •  •  •  ,v0)  .  ■ 

The  converse  of  Lemma  2.2  is  given  by 

Lemma  2.3.  Let  U'mn(z)/V'mn(z )  be  a  Pade  form  of  type  (m,  n)  for  A(z)  €  U,  and  let  D(z), 
D( 0)  4-  0,  be  a  common  divisor  of  U'mn  and  V'mn.  Then  Umn{z)/Vmn{z),  where 

UM(z)  =  U'mn{z)/D(z)  and  V^z)  =  V'^z)/ D(z)  , 

is  also  a  Pade  form  of  type  (m,  n)  for  A(z). 

Proof:  The  arguments  for  the  proof  are  similar  to  those  of  the  proof  of  Lemma  2.2.  a 

The  condition  that  D(0)  ^  0  in  Lemma  2.3  is  necessary  as  can  be  seen  from  the  following 
simple  example. 

Example  2.1:  Let  A(z)  =  1  +  z  +  z2  +  2z3  +  3z4  +  4z3  +  5z6  •  •  a  . 

Clearly,  U'2  A  —  z  and  V'2  A  —  z  —  z1  —  z*  yield  a  Pade  form  of  type  (2,4)  since 

AOO-v'vM  -  (z)  =  o(z7) . 

A  common  divisor  of  U'2ii,(z)  and  V'2,4(z)  is  D(z)  =  z.  However,  the  rational  function  U(z)/V(z), 
where 

U(z)  =  U'2A(z)/D(z)  =  1  and  V(z)  =  V'2A(z)/ D(z)  =  1  -  z  -  z3, 


does  not  yield  a  Pade  form  of  type  (2,4)  since 


. 
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A(z)-V(z)  -  U{z)  =  0(z6).  ■ 


2.2.3.  General  Structure  of  Pade  Forms 

Let  H ^  be  given.  E  det (H^)  =£  0,  then  system  (2.8)  exhibits  a  one-parameter  family  of 
solutions.  Indeed,  if  det(//„^,)  0  and  v0  =  1  in  (2.8),  then  the  solution  of  (2.8)  for  Vmn(z )  and 

of  (2.5)  for  (z)  is  unique.  Pade  forms  for  which  V(0)  =  1,  i.e.,  v0  =  1  are  said  to  be  normal¬ 
ized.  In  addition,  we  have 


Lemma  2.4.  If  det(//m„)  4-  0,  then  the  solutions  U ^  and  V ^  of  (2.5)  and  (2.6)  are  relatively 
prime. 

Proof:  Suppose  the  contrary,  and  let  D(z)  be  a  non-trivial  common  divisor  of  Umn(z)  and  Vm„(z). 

E  £>(0)  =  0,  then  V,„*(0)  =  0,  i.e.,  v0  =  0.  Consequently,  (2.8)  has  only  the  trivial  solution 
and  this  contradicts  the  definition  of  a  Pade  form. 

E  £>(0)  4-  0,  let  £>(z )  =  D(z)/£>(0).  Tlien,  by  Lemma  2.3,  U'(z)  =  U{z)/D(z)  and 
V"(z)  =  V(z)/D(z)  are  solutions  to  systems  (2.5)  and  (2.6).  But  V'(0)  =  V(0),  and  equation  (2.8) 
together  with  the  condition  that  det^^,)  ^  0  then  imply  that  V'(z)  =  V(z),  Thus,  D(z )  must  be  a 
polynomial  of  degree  zero,  which  contradicts  the  initial  assumption  that  D(z)  is  non-trivial.  a 

The  converse  of  Lemma  2.4,  however,  is  not  true,  as  can  be  seen  in  the  following  example. 


Example  2.2:  Consider  again  the  power  series  A(z)  in  example  2.1.  It  can  be  verified  that  a 
Pade  form  of  type  (4,  3)  for  A(z)  is  t/4>3(z)/V4i3(z)  =  (1  -  z  +  z3)  /  (1  -  2z  +  z2),  where 
t/4  3(z)  and  V43(z)  are  relatively  prime.  However,  the  Hankel  matrix 


1  2  3 

2  3  4 

3  4  5 


associated  with  the  Pade  form  of  type  (4,  3)  is  singular,  so  that  det(//43)  =  0.  ■ 
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Given  the  Hankel  matrix  H^,  let  tfm_x  >n_x  be  its  largest  non-singular  principal  submatrix. 
That  is,  X  is  the  smallest  non-negative  integer  for  which  det(//m_x  „_x)  ^0.  From  Lemma  2.4,  a 
Pade  form 


m  -  X 

2>,V 

i-Q 
n-  X 

2v,V 

1-0 


(2.11) 


of  type  (  m-\,n-\  )  for  A(z)  exists,  where  6C,-x,.i-xO)  and  VC-x,«-xO)  are  relatively  prime 
and  vq  ^  0.  The  Pade  form  is  unique  up  to  a  multiplicative  constant  and  satisfies 


and 


'  ’  ’  am-\ 

Vn-X 

“o 

a-n*\  '  °  *  a0 

1 - 

o  • 

—  n  + 1 

-X+ 1 

Vn-X 

V 

-  X 

+n  -  2X 

q 

In  addition, 

AO)  •  VC_x>fl_x0)  “  CC-x,«-xO)  =  zm-x+n~K  +  k+1'2rizi, 


where  tq^O,  if  &<<*>  . 


(2.12) 


(2.13) 


(2.14) 


Theorem  2.5:  Let  X  be  such  that  . _x  is  the  largest  non-singular  principal  submatrix  in 
Hm„  ,  and  let  k  be  determined  by  (2.14).  Then  k  ^  X.  In  addition,  let  1=  min  {  7k, k  } .  Then  a 
basis  for  the  solution  space  of 


• 

' 
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am-n  +  1  •  •  •  +  1 

v*  ‘ 

0 

_ar,  •••  a„+n 

V0 

0 

is  given  by  the  /— X+l  vectors 

{(Cx,  ■  •  '  ,v;,0x)',  •  '  ■  ,  (0,-x,v,‘-x,  •  •  •  .vJ.Ojx.,)'}. 


(2.15) 


Proof:  If  X  =  0,  then  the  theorem  follows  trivially  from  Lemma  2.4.  If  X^O,  then  the  multipli¬ 

cation  of  (2.15)  by 


n  —  X 


X< 


0 


v„_ 


n-k 


o 


0 

1 

Vn-X  •  •  ’  Vq 


yields 


dm  —  «  + 1 


dm  —  X 


0 


■ 

dm  + 1 

V 

■ 

0 

• 

d~m  +n-X 

• 

— 

• 

>■ 

'  1 

■ 

• 

• 

•  •  •  r2X-i-l 

Vo 

0 

(2.16) 


in  the  case  that  k  <  X  ,  or 


’ 
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•  •  •  am 

4- 1 

v„ 

0 

•  •  •  • 

4-n  —  X 

0/-x 

• 

— 

a 

r0 

• 

* 

r0 

r\ 

ro 

r2\-l-\ 

Vo 

0 

(2.17) 


in  the  case  that  k  >  X  .  But  the  largest  non-singular  principal  submatrix  in  (2.16)  has  order  n, 
which  violates  the  assumption  that  X  ^  0  .  Thus,  when  X  ^  0  ,  only  (2.17)  is  possible,  which 
proves  that  k  >  X  . 


From  (2.17),  it  follows  that  the  coeffcient  matrix  of  (2.17)  has  rank  n  +  X  —  /,  and  there¬ 
fore  the  solution  space  of  (2.15)  is  of  dimension  /  —  X  4-  1.  But,  using  (2.13)  and  (2.14),  it  is 
easy  to  see  that  all  the  vectors 

{(Vn-X>  '  "  '  5V0>0\)(  >  '  ’  "  ,  (0/_x,V„_x,  1  •  •  ,Vo,02x-/)<} 


satisfy  (2.17).  ■ 


Corollary  2.6  (Pade):  The  general  solution  of 


and 


«m 

—n 

• 

• 

" 

• 

• 

• 

U0 

On  -n  + 1 

am  + 1 

3m  V, 


'  Qyn  +i 


0 


0 


is  given  by 
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u^z)  =  ^-'(sV)("^V) 

i- 0  i-0 


/-X  n  -  X 


Vm,(z)  =  22^(2a^)(2v;z<), 

i-  O  i-0 


where  l  and  X  are  defined  in  Theorem  2.5  and  <xt ,  0  <  i  <  /— X  ,  is  arbitrary. 

Proof:  Linear  combinations  of  the  basis  vectors  for  the  solution  space  can  be  written  as 


(2.18) 

(2.19) 


Vm.  = 


Vn-X 

■  0/_x  ■ 

• 

v«-x 

«/-x 

; 

+  ••'  +  00 

• 

• 

• 

Vo 

Vo 

°x 

Qzx-z 

Ox 

Oo 

o 

-< 

1 

a 

°x 

a/-x 

02X-i 


v«-x 


Vo 


which  in  polynomial  form  becomes  (2.19).  Furthermore,  using  equation  (2.12)  and  the  results  of 
Theorem  2.5,  it  follows  that 


0/ 

’  am 

- 1 

.0 

i£-k 

• 

Vn-X 

• 

“o 

•  1 

• 

V’O 

o x-i 

a-n  ‘  a0 

T 

1 

for  0  <  i  <  /-X.  Thus,  the  polynomial  form  for  U^z)  corresponding  to  V^(z)  is  given  by 

(2.18).  ■ 


1 
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2.3,  Pade  Fractions 

2.3.1.  Scaled  Pade  Fractions 

Let  for  A(z )  e  U  be  given,  and  let  Hm_y  „_x  be  its  largest  non-singular  principal  subma¬ 
trix.  If  the  rational  function 


'C-\, .-»(*) 


i-0 


2>,V 


i-0 


is  the  Pade  form  of  type  (m-X,n-X)  for  A(z),  then 


(2.20) 


Definition:  The  scaled  Pade  fraction  of  type  (m,n)  fcv  A(z)  is  defined  to  be 

7™(z)  =  ^(z)  /  Tmn(z)  ,  where 

Smn(z)  =  ^  2  (2.21) 

i-0 

L(z)  =  zx2viV.  (2.22) 

i-0 

Theorem  2.7:  The  scaled  Pade  fraction  of  type  (m,n)  for  A(z)  is  a  Pade  form  of  type  ( m,n ) 
for  A(z).  Furthermore,  is  unique  up  to  a  multiplicative  constant. 


Proof:  The  scaled  Pade  fraction  is  obtained  simply  by  setting  a,_x  =£  0  and 

00=01!=  ■  ■  ■  =  o/_x_!  =  0  in  equations  (2.18)  and  (2.19)  of  Corollary  2.6.  ■ 


As  an  immediate  consequence,  we  have 

Corollary  2.8.  det(//mrt)  ^  0  if  and  only  if  rmn(0)  0. 

An  alternative  definition  of  scaled  Pade  fractions  is  given  by 


' 
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Theorem  2.9.  5^(z)/T^(z),  where  Tnn(z )  =£  0  is  the  scaled  Pade  fraction  of  type  ( m,n )  for  A(z) 

if 

(a)  min  {  m  -  d(Smn)  ,  n  -  d(Tmn )  }=  0, 

(b)  GCD(  ,  T™  )  =  for  some  and 

(c)  A(z)  •  7^(z)  -  Smrt(z)  =  0(  z”+"  +  1  ). 

Proof:  The  result  is  an  immediate  consequence  of  Theorem  2.7.  ■ 


Thus,  when  normalized,  the  scaled  Pade  fraction  of  type  (m,n)  for  A(z)  exists  uniquely,  and 
in  addition  satisfies  the  systems  (2.5)  and  (2.6).  This  is  tremendous  advantage  for  our  purposes 
over  other  definitions  of  Pade  fractions  described  below,  since  it  simplifies  the  development  of 
algorithms  in  subsequent  chapters. 

2.3.2.  Non-scaled  Pade  Fractions 

Let  U^iz)  /  Vm*(z)  be  a  Pade  form  of  type  ( m,n )  for  A(z),  and  let 

D(z)  =  CCDiU^  ,Vm). 

Gragg  [GH4 72]  defines  the  Pade  fraction  of  type  (m,n)  for  A(z)  to  be  Pmn(z)  /  Q^iz)  ,  where 

P..(z)  =  Um(z)/D(z) 

=  V„(z)/D(z)  , 

and  in  addition  Qmn  (0)  =  1. 


P^(z)  =  z  x  Smn(z)  =  lfn ,_XiW_x(z) 


That  is, 


■ 
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Qmn  00  =  T^iz)  =  VC_Xi(t_x(z)  , 

where  I/’ _x  „_x(z)  and  VC_X  „_x(z)  are  defined  in  (2.20).  Thus,  the  Pade  fraction  of  type  ( m,n ) 
always  exists  and  is  unique.  However,  /^(z)  /  (z)  is  a  member  of  the  family  of  Pade  forms  of 

type  (m,n)  for  A(z)  given  by  (2.18)  and  (2.19)  only  if  /  =  2\  ,  or  equivalently,  if  k  >  2\  in  equa¬ 
tion  (2.14).  Thus,  in  the  case  that  k  >  2A  ,  the  Pade  fraction  is  obtained  from  the  Pade  form  by 
setting  a0  =  l,ax  =  •  •  •  =  at_k  =  0  in  equation  (2.18)  and  (2.19).  In  the  remaining  case  that 
k  <  2\  ,  the  Pade  fraction  Pmn(z )  /  0^(z)  does  not  satisfy  the  systems  (2.5)  and  (2.6). 

Baker’s  [BAK15]  perspective  is  different  in  a  very  subtle  way.  He  requires  that 
Pmn(z)  /  Qmn 00  satisfies  (2.5)  and  (2.6)  and  in  addition  0^(0)  =  1  and  GCD(Pmn,QmA )  =  1  .  But 
once  again,  from  Theorem  2.5,  this  is  possible  only  if  k  >  2\  ,  in  which  case  Pmn(z)  /  Qmn{z )  is 
obtained  by  setting  -  1,0^  =  •  •  ■  =  a,_x  =  0  in  the  equations  (2.18)  and  (2.19).  Thus, 
from  Baker’s  point  of  view,  a  Pade  fraction  may  not  exist,  but  whenever  it  does,  it  is  unique  and 
satisfies  systems  (2.5)  and  (2.6). 


. 


CHAPTER  3 


COMPUTATION  OF  OFF -DIAGONAL  SCALED  PADE  FRACTIONS 


3.1.  Introduction 

In  this  chapter,  derived  is  a  new  relationship  between  three  scaled  Pade  fractions  lying  along 
an  off-diagonal  path  of  the  Pade  table  for  a  power  series  A(z).  This  relationship  is  described  in 
Section  3.2  and  is  used  in  Section  3.3  to  develop  an  algorithm  to  iteratively  compute  a  sequence  of 
scaled  Pade  fractions  along  the  off-diagonal  path.  However,  the  algorithm  is  not  totally  iterative, 
since  the  relationship  recursively  involves  a  scaled  Pade  fraction  of  a  different  power  series  A(z), 
computed  from  A  (z).  By  doubling  the  step-size  along  the  off-diagonal  path  at  each  iteration,  it  is 
shown  that  the  complexity  of  the  algorithm  for  computing  the  n-th  scaled  Pade  fraction  along  the 
path  is  (n  log2n),  assuming  fast  polynomial  methods  are  used. 

The  algorithm  of  Section  3.3  can  be  used  to  compute  scaled  Pade  fractions  for  the  quotient 
-A(z)/£(z)  of  two  power  series  A(z)  and  fl(z).  This  can  be  accomplished  by  formally  computing 
the  inverse  of  B(z),  multiplying  it  by  -A(z),  and  then  applying  the  algorithm  to  the  result.  The 
dominating  cost  of  this  procedure  is  the  cost  of  computing  the  inverse  of  Z?(z).  This  cost,  however, 
can  be  eliminated  by  an  obvious  modification  of  the  algorithm,  which  is  given  in  Section  3.4. 

If  the  step-size  along  the  off-diagonal  path  is  shortened,  all  scaled  Pade  fractions  along  the 
path  can  be  computed.  The  recursive  call  alluded  to  earlier  then  becomes  trivial,  and  the  resulting 
algorithm,  given  in  Section  3.5,  becomes  truly  iterative.  In  the  case  that  the  path  is  along  the  diag¬ 
onal,  the  algorithm  is  identical  to  the  one  given  by  Cabay  and  Kao  [CAB 83].  This  algorithm  can¬ 
not  take  advantage  of  fast  polynomial  methods,  and  has  complexity  0(n2). 
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3.2.  Preliminary  Results 

The  scaled  Pade  fractions  can  be  arranged  in  a  doubly  infinite  array  as  follows: 
Definition:  The  collection  of  all  scaled  Pade  fractions  of  type  (m,n)  for  A(z)  e  U,  given  by 


r  (a) 


7-i,-i 

7-1,0 

7-1,1 

•  •  ■  7-i,* 

7o,-i 

7o,o 

7o,i 

•  ■  ■  7o^i  • 

7i.-i 

7i,o 

7i,i 

•  •  •  7l,n  • 

7m, -1 

7m  ,0 

7m  ,1 

•  •  •  7m,*  • 

(3.1) 


is  called  the  (extended)  scaled  Pade  table1  for  A(z).  The  modifier  extended  denotes  the  inclusion 
of  the  first  row  and  column  in  the  table  which  are  defined  as  follows: 

Sm -i(z)  =  -  f  and  Tm_x{z)  =  0  ,  for  n>  -1,  and 
S_ln(z)  =  0  and  7_ln(z)  —  —  z*  ,  for  0.  ■ 

For  computational  purposes,  define  the  ^/-truncated,  scaled  Pade  table  for  A(z)  to  be 


r*(A) 


T-i ,-i  *y-i,o  7-i,i  •  •  •  7-i w 

70, -i  7o,o  7o,i  '  •  -  7ow 

71,  — l  7i,o  7i  ,i  •  •  ’  71,* 


(3.2) 


Let  m  and  n  be  non-negative  integers  such  that  n  >  N  .  The  next  few  results  are  concerned 


!Note  that  the  (extended)  scaled  Pade  table  corresponds  to  the  usual  definition  of  the  Pade 
table  given,  for  example,  in  Gragg  [GRA12\  with  the  exception  that  the  scaling  factor  z*  in  (2.21) 
and  (2.22)  does  not  appear. 


' 
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with  the  construction  of  the  scaled  Pade  fraction  -y^z),  given  that  ^(A)  already  exists. 
Without  loss  of  generality,  assume  that  m  ^  n  (otherwise,  the  same  arguments  can  be  applied  to 
1  /A(z)  eU).  Let 

M  =  N  +  (m  -  n)  .  (3.3) 

Then,  (m,  n )  and  ( M ,  N)  both  lie  along  the  (m-n)th  off-diagonal  path  of  the  scaled  Pade  table 
T(A)  (see  Figure  3.1  ). 


N 


For  the  scaled  Pade  fraction  of  type  (M ,N)  for  A(z),  equation  (2.3)  becomes 

+N  +  *w+l 


A(z)  ■  Tmn(z)  Smn(z)  —  z 


*l(z)  , 


(3.4) 


I  \ 


1 1 
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where  jxWA,  >  0,  and  /?i(0)  ^  0  if  m,*  <  00  .  Let 

Zx*"=  GCD(Sm,Tu„)  .  (3.5) 

To  construct  "ym„(z)  ,  two  separate  cases,  \lmn  >  ( n-N )  and  <  ( n-N ),  arise.  These  two 
cases  are  considered  seperately  in  Theorem  3.1  and  Theorem  3.4. 

Theorem  3.1:  If  y^MN  >  (n-N)  in  (3.4),  then  the  scaled  Pade  fraction 

=  5™ (z)/T„„(z)  is  given  by 

S„„(z)  =  2*-*  Sm(z)  (3.6) 

Tmn(z)  =  ?-*  Tun(z)  .  (3.7) 

Proof:  Qearly  from  relation  (3.3),  it  follows  that 

^(Smji)  =  n—N+d(SUft f) 

^  n +M—N 
=  m. 

Similarly,  d(Tm„)  <  n  .  Furthermore,  it  is  clear  that  min{  m-d(Smn),  n-bij^  }  =  0  from  the 
equality,  min{  M - d(SMN) ,M - d(TMN)  }  =  0.  Using  the  fact  that  (n-N)  <  \iMN,  it  follows  that 

A(z)Tmn(z)  -  S^z)  =  zn~N  (A(z)Tmn(z)  -  SMN(z)) 

=  (Z""''*'M'"'«1(Z)) 

=  z"*"x,‘*K*1fi1(Z) 

=  . 

Finally,  GC£>(S„,  ,Tm.)  =  CCD(Sm,Tm)  =  . 

Consequently,  y mn(z)  -  Smn(z)  /  Tmn(z)  given  by  (3.6)  and  (3.7)  is  the  scaled  Pade  fraction 


of  type  (m,n)  for  A(z)  . 
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For  the  case  [iMN  <  ( n—N )  ,  let 

A/*  —  M  —  \MN  —  1  ,  (3.8) 

N*  =  N  -  \MN  —  1  (3.9) 

(see  Figure  3.1  ).  dearly,  the  scaled  fraction  7wv(z)  for  Mz)  satisfies 

A00  •  WOO  ~  =  zM^N^R0(z)  ,  (3.10) 

where  i?0(0)  ^  0  .  In  addition,  by  Theorem  2.7  and  by  the  definition  of  M *  and  N*,  it  follows  that 

Smn(z)  /  Tmn(z)  SM*N*(z)  /  Tm*n»(z )  .  (3.11) 

From  the  two  unit  power  series,  i?i(z)  and  i?0(z)  given  by  (3.4)  and  (3.10),  respectively, 
construct  the  unique  power  series 

A(z)=/?0(z)/JR1(z).  (3.12) 

Associated  with  A  (z)  eU ,  let 

m  =  n  -  N  +  \MN  ,  (3.13) 

n  =  n  —  N  —  \lmn  —  1  .  (3.14) 

Now  assume  that  the  n-truncated  scaled  Pade  table,  T-(A),  for  A(z)  has  been  constructed,  as  well. 
That  is,  assume  that  the  ( m,n)  scaled  Pade  fraction,  y--(z)  =  S--(z )  /  T--(z),  such  that 

A(z)  ■  T-(z)  -  S--(z)  =  0(  z’*"-*1  )  (3.15) 

is  available. 

The  scaled  Pade  fractions  y MN{z. ),  *yM  vOO  'Y«;r(z)  provide  suffiaent  information  to 
obtain  directly  a  Pade  form  of  type  (m,n)  for  A(z)  .  This  Pade  form  is  constructed  in  Lemma  3.2 
below.  It  is  shown,  later  in  Theorem  3.4,  that  this  form  is  also  the  scaled  Pade  fraction  of  type 
( m,n )  for  A(z). 


% 
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Lemma  3.2:  Let  \iMN  <  ( n—N )  in  equation  (3.4),  and  let 

S„„(z)  =  z-lulSm(z)  S-,(z)  -  z“  S„.M.(z)  f--(z)  (3.16) 

T™(z)  =  z)  S-;(z)  -  z“  r„v(z)  f--(z),  (3.17) 

where  a  =  \MN  +  \iMN  +  2  .  Then  Smn(z)  /  T^z)  is  a  Pade  form  of  type  ( m,n )  for  A(z). 

Proof: 

a(S„)  =  max{a(z_l*"S„„S-;),  3(z“SJ(Vfs;)} 

<  max{  —  Xw/Vr  +  M  +  m,  a  +  M*  +  F} 

=  rnax{  —\MN  +  M  4-  (n  —  N  4-  XWN), 

(^MN  +  \^MN  +  2)  +  (M  —  \MN  —  1)  +  (n  —  N  —  H'MN  —  1)} 

=  max{m,  m) 

=  m. 

Similarly,  a(7mn)  <  n.  Moreover,  the  fact  that  SMN(z)  /  Tmn(z),  Sm.n.(z)  /  Tm.n*(z)  and 
S--(z)  /  f--(z )  are  tiie  scaled  Pade  fractions  yields  immediately  that 

min{rn  -  a(Smfl)  ,  n  -  a(rmil)}  =  0. 

Furthermore,  from  (3.4),  (3.10),  (3.12)  and  (3.15),  it  follows  that 

A(z)  •  Tmn(z)  -  Smn(z)  =  A{z)  {z~^Tun{z)S--{z)  -  z“  T„v(z)f--(z)} 

-  {z~ KmSMN (z)S--(z)  -  f  5wV.(z)f--(z)} 

=  ^"^--(z)  {A(z)r^(z)  -  SM„(z)} 

-  yoMWr.vW  -  suVW) 

=  z""'Si;„-(z){Ri(z)  -  Z“  f^zJWfz)  z"'*"%1} 
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_  zM+N  +  *lM^-KW+1 


Sxk)  -  Ro(z)  T--{z)} 


=  -Rx{z)  zm+n+^km+1 


{  Tmn(Z)  Ro(Z)  /  *l(z)  "  S--(Z)  } 


=  -Ri(z)  z 


,M+N+tlAfr 


0(zm+H  +  1 ) 


Q  (  ZM  + N  + ^  ATT  KAfA  1 +  m  +  "  + 1  j 


=  0(zm  +  n+1).  B 

In  order  that  Smn  /  in  (3.16)  and  (3.17)  be  a  scaled  Pade  fraction  of  type  (m,  n ),  it 
remains  to  show  that  GCD(Smn,Tmn)  =  zm  for  some  .  With  this  intent,  consider  again  the 
(in,  n)-th  entry  of  F-(A),  and  let 


m  =  m  —  X--  —  1 

m  n 

n  =m  -  X--  -  1, 

where 

zX”"  =  GCD(S--,  f--). 

^  m n  7  win  / 

The  scaled  Pade  fraction  y_»  _*(z)  for  A(z)  satisfies 

AT -  5-.-.  =  0£(z"%"*+1) 

m  n  m  n  ' 

_  x  m  +  n  -  2A~“  -1.  * 

=  <9*(z  mfl  ),2 


and,  in  addition, 

5--(z)  S-*~*(z) 

mn  V  /  m  n  v  7 

f--(z)  f_._.(z)  ‘ 

mnv  /  m  n  '  ' 


(3.18) 

(3.19) 


(3.20) 


(3.21) 


(3.22) 


2  R  (z)  =  0E  ( zk  )  means  that  R  (z)  is  a  power  series  whose  first  non-zero  coefficient  is  the 
coefficient  of  zk,  exactly. 


. 
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Lemma  3.3:  Let 


W*)  4-  .-*(0  -  •  „•«  r_.  _.(z) 

(3.23) 

=  z'‘*"'T„„(z)4.-.(z)  -  z*  r„.„.(z)f_.;.(z) 

(3.24) 

where 

«  =  Vw  +  +  2  .  Then  -ym.  ^.(z)  =  $„•„•(*) /r.v(z)  is  a  Pade 

form  of  type 

(m  ,  n 

)  for  A(z),  where 

m*  =  m  —  X-  -  —  1 

m  n 

(3.25) 

n  =  n  -  X---  1. 

(3.26) 

Proof:  The  proof  is  identical  to  the  proof  of  Lemma  3.2.  ■ 


Theorem  3.4:  Let  | N  <  (n  -  N)  in  equation  (3.4),  and  let 

=  Smniz)  ' 

be  defined  by  (3.16)  and  (3.17).  Then  y mn(z)  is  a  scaled  Pade  fraction  of  type  (m,  n )  for 
A(z). 

Proof:  Let 

Gmn(z)  -  GCD{Smn,Tmn). 

We  first  show  that  ^(G^)  ^  X--,  where  X--  is  given  by  (3.20).  Suppose  that 
3(Gm*)  >  X-  -  ,  and  proceed  by  contradiction.  Let 

um.  „.(z)  =  z'<°”"’'k™~1  ^(Z)  /  G„(z)  , 

V..  „.(z)  =  T„(i)  /  G„(z), 

where  m*  and  n*  are  given  by  (3.25)  and  (3.26).  Then  d(Um*n •)  ^  m* ,  3(Vm»n*)  <  n* ,  and 

*(*)  •  v.-m  -  <v.-w 


' 


' 
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=  zs(c”")“‘""_1{/4(2)r.„(z)  -  S.,(z)}/G.,(z) 

_  0(z(_X"r"  _1)  +  (m  +"  +  1)) 

,  m *+«*  +  X  ~~+ 1 

=  0(z  mn  ). 

Thus,  Um'n>(z)  /  Vm*n>(z)  is  a  Pade  form  of  type  {m*  ,n)  for  A(z). 

But  Sm»n*(z)  /  Tm*n»(z),  given  by  Lemma  3.3,  is  also  a  Pade  form  of  type  ( m  ,  n  ).  Then, 
from  the  general  structure  of  Pade  forms  given  in  equations  (2.18)  and  (2.19)  of  Corollary  2.6,  it 
follows  that 


S  *  »(z)  /  T  •  »(z) 

in  n  '  '  in  m  '  ' 


=  U  .  .(z)  /  V  .  .(z) 
Smn(z)  /  Tnn{z), 


or,  equivalently,  that 


S*  *(z)  Tnn(z)  -  Sm „ (z)  T  •  *(z)  =  0. 


(3.27) 


Replacing  Smn(z)  /  Tm„(z)  and  in  equation  (3.27)  by  the  expanded  forms  (3.16), 

(3.17),  (3.23)  and  (3.24),  it  follows  that 


z“rx*"'{SM„(z)  T„  v.(z)  -  S„ .„.(z)  Wz)HS„-„-{z)  T'-.-.fz)  -  5„.„.(z)  T-. _.(z)}  =  0  . 


Thus,  either  S„ K(z)  /  Tu „(z)  =  Su->l.(z)  /  T^.^-iz)  ,  or  S--(z)  /  T-^z)  =  S_»  _.(z)  /  r_.  _.(z)  , 
which  contradicts  (3.11)  and  (3.22). 


Thus,  a(G„.)  s  X-- .  But, 

z*™"  =  GCD(S-  -  ,  T- 

which  implies  that  zXmn  divides  both  Sm„(z)  and  Tmn(z )  in  equations  (3.16)  and  (3.17).  That  is,  zXmn 
divides  Gm  „  (z)  ,  and  consequently 


J 
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Gmn(z)  =  zX""  . 

Thus,  Smn(z)  /  Tmn(z)  given  in  Lemma  3.2  is  not  only  a  Pade  form  of  type  (m,n)  for  A(z),  but  also 
the  scaled  Pade  fraction  of  type  (m,n)  for  A(z),  where 

GCD(S M  ,  Tmn )  =  zX- 

and  .  ■ 

Theorem  3.5:  Let  <  ( n-N)  in  (3.4)  and  let  7m„(z)  =  Smn(z)  /Tmfl(z)  be  the  scaled  Pade 
fraction  of  type  ( m,n )  for  A(z) .  If 

m  =  m  -  Am*  -  1  (3.28) 

and 

n  =  n  —  \mn  —  1  ,  (3.29) 

where 

=  GCD(S„.  ,  r„„)  ,  (3.30) 

then  the  scaled  Pade  fraction  of  type  (m*,n*)  for  A(z)  is  =  1  v(7)>  wiiere 

5m*fl»(z)  and  Tm»n*(z)  are  given  by  equation  (3.23)  and  (3.24). 

Proof:  The  theorem  follows  using  arguments  identical  to  those  of  proof  of  Theorem  3.4  ,  and 

using  the  results  of  Lemma  3.3.  ■ 

A  simple  example  for  the  off-diagonal  computation  is  presented. 

Example:  Let  A(z)  =  1  +  z4  +  z5  +  z9  +  z10  +  2z1!  +  •  ■  •  .  This  example  constructs  the 

scaled  Pade  fraction  7?i6 (z)  =  S1(>(z)  /  T7  6(z)  of  type  (7,6)  for  A(z).  Since  m  -  n  =  1,  the  con¬ 
struction  proceeds  along  the  1st  off-diagonal  path  of  the  scaled  Pade  table  T(A) . 

Assume  that  T3(A)  is  already  available,  from  which  it  can  be  determined  that  the  scaled 
Pade  fraction  -y4 _3(z)  of  type  (M ,N)  =  (4,3)  is  given  by 
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■^4,3  (0  _  1  ~  Z  +  Z7  ~  2?  +  Z* 
^4.3  (^)  1  ~  Z  +  Z2  ~  Z3 


From  (3.4),  the  residual  for  S43(z)  /  r43(z)  is  given  by 

A(z)  74>3(z)  -  54>3(z)  =  ^(z)  z4+3  +  1  , 
where  R\(z)  =  —  1  +  z  —  z5  +  2z7  +  •  •  •  . 


(3.31) 


(3.32) 


Consequently,  equation  (3.32)  yields  that  jx4  3  =  0.  Since  p,4>3  <  n  —  N ,  Theorem  3.4  is  applica¬ 
ble.  Observe  that 


zMS  =  ?  =  gcd(s<i3  ,  r4>3), 

and  consequently  the  predecessor  of  *y4>3(z)  along  the  1st  off-diagonal  path  is  yM*N*(z)  =  73,2(2). 
Therefore,  73  2(z)  is  contained  in  T3(A),  and  is  found  to  be 


73,200  = 


^3,2^) 

Tv(z) 


z2  ’ 


The  residual  for  S3  2(z)  /  r3>2(z)  is  given  from 

A(z)  ^3,2(0  “  S3i2(z)  =  /?0(z)  z3+2  +  1  ,  (3.33) 

where  /?0(z)  =  1  +  z  +  z5  +  z6  +  •  •  . 

The  two  residuals  /?i(z)  and  R0(z)  give  the  residual  power  series  A  (z)e  U,  where 
A(z)  =  tf0(r)  /  i?i(z) 

=  -1  -  2z  -  2z2  -  2z3  -  •  •  •  .  (3.34) 


Using  equation  (3.13)  and  (3.14), 

in  —  n  —  N  +  \MN  =  3 
n  =  n-  N-\LMN  -1  =  2, 

and  in  order  to  apply  Theorem  3.4,  it  is  therefore  required  to  obtain  the  scaled  Pade  fraction 
73  2(z)  of  type  (3,2)  and  its  predecessor  for  A(z).  Assuming  that  T3(A)  is  available,  from  it  can  be 
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determined  that 


-  ,  \  $3.2(2) 

73,2(2)  =  y— 
*3,2(2) 


— z  —  z 


z 2 


(3.35) 


Similarly,  the  predecessor  of  73,2(2)  along  the  first  off-diagonal  path  of  T3(A)  is  given  by 
71,0(2)  =  5li0(z)  /  f1)0(z)  =  (-1  -  2z)/l. 


Now  by  applying  the  formulae  (3.16)  and  (3.17),  the  (7,6)  entry  of  T(A)  is  computed  by 


$7,6(2)  =  (1  -  z  +  z2  -  z3  +  z4)-(-z  -  z2)  -  z*  -  z2)  , 

77>6(z)  =  (1  -  z  +  z2  -  z3)-(-z  -  z2)  -  z*  (z2)-(z  -  z2)  , 


where  a  —  X4>3  +  p,4,3  +  2  =  2.  Thus, 

77,6(2)  =  $7,6(2)  /  T7fi(z)  =  z(-l  -  z4)  / z(  1  +  z!).  (3.36) 

Note  that  71,0(2)  is  not  used  for  computing  77>6(2).  It  is  used  instead  in  formulae  (3.23)  and  (3.24) 
for  computing  the  predecessor  of  7?  6(z) .  Since 

=  2*’”  =  z1  =  GCD(Sxl  ,  f„), 

it  is  known  that  the  predecessor  is  7j,4(z),  which  is  determined  by  (3.23)  and  (3.24)  to  be 

y  (z\  =  554^z)  =  C"1 *...+  2"  -  z3  -  2 z5)  B 
5,4  ^5,4(2)  (— 1  —  z  +  z2  —  z3  +  z4) 


3.3.  Fast  Off-diagonal  Algorithm  for  a  Single  Power  Series 
3.3.1.  The  Algorithm 

The  algorithm  given  in  this  section  constructs  the  scaled  Pade  fraction  ymn  of  type  ( m,n )  for 
A(z)  in  a  quadratic  fashion.  The  iteration  assumes  the  existence  of 

7a !n(z)  ~  $mn(z)  !  Tmn(z)  ,  (3.37) 

where  M  =  N+(m  —  n),  and  where  for  notational  convenience  we  set 
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•Si  —  Smn(z) 


and 


^1  —  Tmn(z)  . 


Also  assumed  to  exist  is  the  predecessor 

V.'W  =  WW/!i,v(!) 

of  7 MN (z)  on  the  the  (m-/i)*th  off-diagonal  path  of  T(A),  where 

If  1  , 

AT  =  N  -  XM*  -  1 

and 

zx*"  =  GCD(Sun,Tu„)  . 


Again  for  notational  convenience,  we  set 


and 


So  = 


S  .  .(z) 

A/  A1 


To  advance  the  solution  from  N  to  N  +  s  (that  is,  to  construct  7^ +  d,v+j(z));  where 
step  size,  the  algorithm  first  computes  \iMN  such  that 


A(z)T \  -  Sx  mod  “+N+2'+^'  =  /  +  »  +  ^+1  /jl(2)  , 


where 


zXm  =  GCD(  Sx  ,  7\  ) 


and  fli(0)  *  o  if  <  2s  + 


(3.38) 

(3.39) 

(3.40) 


(3.41) 

(3.42) 
s  is  the 


If  \lmn  >  s  ,  then  7W+J,N  +  J(r)  is  constructed  trivially  by  means  of  Theorem  3.1.  Otherwise, 


1 


Theorems  3.4  and  3.5  are  applied. 


ALGORITHM  1:  OFFDIAG 

INPUT:  A,  m,  n  ,  where 

(1)  m  and  n  are  non-negative  integers  with  m  >  n  ,  and 

(2)  A  is  a  unit  power  series  .  (  Note  that  only  A  mod  zm+"+1  is  required  ). 


OUTPUT: 


where 


(1)  /  Tx  is  the  scaled  Pade  fraction  of  type  (m,  n )  for  A,  and 

(2)  S0  /  T0  is  the  scaled  Pade  fraction  of  type  (m-Xrin-l,  for  A,  given  that 

zkm  =  GCD(S1,Tl)  . 

Step  1:  #  Initialization  # 

i  -  —  1 
M  -  (m  —  n) 
jV  -  0 


Si  S0 

A(z)  modzM  +  1  -z"-1  ' 

T'i  To 

1  0 

Step  2:  #  Calculation  of  step-size  # 


i  «-  i  A  1 

s  -  min  {2?  -  N  ,  n  —  N} 


Step  3: 


#  Termination  criterion  # 
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If  5  =  0  then  exit 

Step  4:  #  Calculation  of  scaling  factor  for  Sx  /  Tx  # 

Determine  \MN  such  that  zXmv  =  GCD{  SX,TX  ) 

Step  5:  #  Computation  of  residual  for  yMN  =  Sx  /  Tx  # 

Compute  and  such  that 

(ATX-SX)  mod  zM  +  N+2,+  K^1  =  +  Rif 

where  Rx  (0)  ^  0  if  \iMN  <  7s  +  \MN . 

Step  6:  #  Identification  of  Cases  # 

if  ^MN  —  s 

then  #  Case  of  Theorem  3.1  # 


Sx  Sq 

S  i  So 

V  o' 

i _ 

A  r0 

0  1 

go  to  step  12 

else  #  Case  of  Theorem  3.4  and  3.5  # 
go  to  step  7 

Step  7:  #  Calculation  of  degrees  for  residual  scaled  Pade  fractions  # 

m  -  s  +  \MN 
n  -  5  —  |xMjV  —  1 

Step  8:  #  Computation  of  residual  for  -  S0  /  T0  # 


Compute  R0  such  that 


' 
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(A(z)To  So) 


.  M  +  N  +  m+  n  —  2X  »*’ 

mod  z  = 


R 


O  5 


where  /?0(0)  =£  0. 


Step  9:  #  Computation  of  residual  power  series  # 

A(z)  -  i?o(2)  /  /?i(z)  (w0<7  z«  +  «  + 

Step  10:  #  Computation  of  residual  scaled  Pade  fractions  # 


Si  S0 
f i  f0 


-  OFFDIAG  ( A(z),  m,  n  ) 


Step  11:  #  Advancement  of  scaled  Pade  fraction  computation  # 


s ! 

5o‘ 

5i 

So 

7i 

T-i 

T0 

0 

_zX^'+tiAf/+ 


Si 

50 

fi 

fo 

. 

Step  12:  #  Calculation  of  degrees  of  Si  /  Ji  # 


N  -  N  +  s 
M  -  M  +  s 
go  to  step  2  e 


3.3.2.  Proof  of  algorithm  validity 

Let  yi'V(z)  be  the  scaled  Pade  fraction  of  type  (m,n)  for  A(z),  where  m  >  n  >  —1  .  As  in 
(3.2),  define 

I*(A)  -  {  y^(z)  |m>n,S>n>-l}  (3.43) 

to  be  the  N-truncated  Pade  table  for  A(z),  with  the  additional  provision  that  m  >  n  . 
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Also,  define 


=  {I\(A)  for  all  AeU} 


(3.44) 


for  N  =  0,1, 

The  proof  proceeds  by  induction  on  for  i  =  0,1,2,  •  ■  .  Qearly  the  algorithm  is  correct 
for  Oq. 

It  is  shown  that  if  the  algorithm  is  correct  for  ,  then  it  is  also  correct  for  Ct^+i ,  where 


Let  A (z)  e  U  be  an  arbitrary  unit  power  series.  It  is  shown  that  the  algorithm  correctly  computes 
,  where  m  >  n  and  —  1  <  n  <  N+L  ,  given  that  it  is  correct  for  (1N  . 

If  —  1  <  n  ^  N  ,  then  *y^(z)e  FN(A)e  ClN  ,  and  by  the  inductive  hypothesis  the  algorithm 
computes  it  correctly. 

If  N  <  n  <  N  +  L  ,  let  M  -  N  +  (m-n).  By  the  inductive  hypothesis,  the  algorithm  (using 
logjiV  iterations)  correctly  computes  yify(z)  =  Sj£j(z) 1  Fj£)(z)  e  rw(A)  e  ClN  ,  and  its  predecessor 
TJw'*)Jv*(z)  •  ^  remains  to  show  that  one  more  iteration  correctly  gives  yj^z)  and  its  predecessor 


Let  s  =  n  -  N  be  given  by  step  2  of  the  algorithm  and  let  z =  GCD(S[$  ,  7^)  be  given 
by  step  4.  Then  step  5  yields  |xWN  >  0  and  a  polynomial  R\(z)  such  that 


(  A(z)Tj$(z)  -  S^X1)  )  ™>d  z 


M  +  N  +  2j  +  X  1 


^Rl(z)zM+N+^+1 


(3.46) 


where  ^  0  ,  and  f?i(0)  ^  0  if  iiMN  <  2s  +  \MN  .  Then  two  cases  arise. 


' 


■ 
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If  P- mn  —  s  ,  then  by  Theorem  3.1,  step  5  correctly  yields  *y£V(z)  .  Furthermore,  from  equa¬ 
tions  (3.6)  and  (3.7), 

zx™  =  gcd(s!$(z)  ,  7<i>(j)) 

_  ZKW+  1 

Consequently,  the  predecessor  of  y^(z)  is  given  by 

■  (3.47) 

since 


and 


B*=»-X„-l=Ar  -  w  -  l  =  n"  . 


If  I^mn  <  z  .  then  step  7  gives 

n  =  s  -  [iMN  -  1  <  N 


and 


m  =  s  +  \MN  >  n  . 

The  polynomial  i?:(z),  computed  in  step  5,  therefore  satisfies  /^(O)  ^  0  and  is  of  degree  m  4-  n  . 
In  addition,  step  8  produces  a  polynomial  i?0(z)  of  degree  m  4-  n  such  that  /?0(0)  0  .  Conse¬ 

quently,  enough  terms  are  available  in  Z?i(z)  and  tf0(z)  to  compute,  in  step  9,  the  residual  unit 
power  series  A(z)  mod  z"  +  "  +  1  defined  by  (3.12).  By  the  inductive  hypothesis,  step  10  there¬ 
fore  correctly  computes  y^£(z)  e  rw(A)  t  ClN  and  its  predecessor.  It  now  follows  by  Theorems 
3.4  and  3.5  that  step  11  correctly  computes  yi'V(z)  and  its  predecessor.  ■ 


,  4 
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3.3.3.  Cost  Analysis 

Let  C(m,n)  be  the  cost  of  computing  the  scaled  Pade  fraction  of  type  (m,n)  and  its  predeces¬ 
sor  for  an  arbitrary  power  series  A(z)  e  U  using  Algorithm  1.  For  the  sake  of  simplicity,  assume 
0  <  n  <  m  <2n .  The  case  that  m  >  2n  is  considered  later.  In  this  section,  asymptotic  esti¬ 
mates  of  C(m,n)  are  derived  by  counting  the  number  of  operations  (additions,  subtractions,  multi¬ 
plications  and  divisions  in  F)  performed  by  the  algorithm,  A  detailed  cost  analysis  and  an  imple¬ 
mentation  of  the  algorithm  are  described  by  Verheijen  [VEP83],  who  compares  Algorithm  1  with 
other  algorithms  for  calculating  Pade  fractions. 

When  obtaining  the  asymptotic  cost  estimates,  it  is  assumed  that  the  algorithm  makes  use  of 
fast  methods  for  polynomial  and  power  series  arithmetic.  Using  fast  Fourier  transforms,  two  poly¬ 
nomials  with  coefficients  in  the  field  F  and  of  degree  M  and  N,  respectively,  can  be  multiplicated 
in  0((M  +  N)  log (M  +  N))  operations  in  F.  Using  Newton’s  method  and  fast  multiplication  of 
polynomials,  the  first  N  terms  of  the  quotient  of  two  power  series  in  U  can  be  obtained  in 
0(N\ogN)  operations  in  F.  These  and  other  fast  methods  for  polynomial  and  power  series 
arithmetic  are  described,  for  example,  in  Aho,  Ho  per  aft  and  Ullman  [AH  01  A]  and  in  Lips  on 
[ZiP  81]. 


Let  k  =  [log  n] .  Then  it  is  easy  to  verify  that  the  algorithm  terminates  after  k  iterations, 
and  that  after  the  execution  of  step  2  of  iteration  i 


s 


1 

2'"1 

n  —  2k~1 


,  i  =  0 
,  0  <  i  <  k 
,  i  —  k. 


(3.48) 


Consequently,  during  iteration  i, 


fO  ,/  =  0 

12'-1  ,  otherwise, 


(3.49) 


and 


■ 
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M  =  N  +  (m-n)  .  (3.50) 

Assume  that  during  iteration  i,  the  scaled  Pade  fraction  Sx/  Tx  of  type  (M,N)  and  its  predecessor 
SQ  /  Tq  are  available. 

The  first  nontrivial  step  requires  the  computation  of  the  residual  ^(z)  in  step  5.  Let 

M 

*\(?)  =  (3.51) 

J- o 


and 


AT+2i-l 

A2(z)=  2  aM+]+ •  (3.52) 

j-o 

With  \MN  given  by  step  4,  it  is  known  that 

A(z)  Tx  -  Sx  mod  zM+N+*+^ 

=  Ax(z)  Tx  +  zm^1A2(z)  Tx  -  Sx  mod  zM+N+2s+x^1 
=  0(  zM+N  +  1).  (3.53) 


Since  AXTX  and  Sx  are  both  of  degree  at  most  M  +  N,  then  and  /?i(z)  can  be 

obtained  directly  from  zM  +  lA2Tx.  The  product  A2TX  is  a  polynomial  of  at  most  degree 
2N  +  2s  -  1  <  AN .  Thus,  using  fast  polynomial  multiplication,  step  5  can  be  executed  in 
O  (N  log  N)  operations. 

If  it  is  determined  as  a  result  of  step  5  that  \iMN  >  s ,  then  step  6  is  performed  trivially  and 
the  iteration  is  complete.  Otherwise,  the  algorithm  continues  in  step  8  with  the  calculation  of  the 
residual  R0  of  the  predecessor  scaled  Pade  fraction  S0  /  T0.  Making  the  same  observations  as  in 
step  4,  it  follows  that  R0  can  be  obtained  from  the  product  of  A3  70,  where 

N+m+n-X^ 

A  3  —  X  au (3.54) 

j-0 

and  T0  is  a  polynomial  of  degree  at  most  N  —  \MN  —  1 .  Since  the  degree  of  the  product  is 
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bounded  by  2N  +  m  +  n  —  2XMN  —  1  <  4V,  step  8  can  be  executed  in  0(N  log  N )  operations. 

Hie  computation  of  the  residual  power  series  A  (z)  mod  z™  +  "  +  1  in  step  9  requires  the 
computation  of  the  first  in  +  n  =  2s  +  \MN  —  —1^3N  terms  of  the  quotient  of  R0  /  Rx. 

Again,  this  may  be  performed  in  0(N  log  N)  operations. 

In  step  8,  the  recursive  call  of  Algorithm  1  in  order  to  compute  the  scaled  Pade  fraction  of 
type  (in,  if)  for  A(z),  requires  C(m,  if)  operations  by  assumption.  For  later  purposes,  it  is 
important  to  observe  that  0  <  n  N  and  that  n  <  in  <  2 N. 

The  final  non-trivial  step  requires  eight  polynomial  multiplications  to  obtain  the  scaled  Pade 
fraction  of  type  (M  +  s,  N+s)  and  its  predecessor  for  A(z).  Snce  M  ^  2N,  each  of  the  poly¬ 
nomial  products  are  of  degree  at  most  M  +  s  <  3N.  Consequently,  step  8  an  be  executed  in 
0(N  log  N)  operations. 

It  is  an  easy  matter  to  show  that 

C(mu  rh)  <  C(m2,  n2),  (3.55) 

whenever  mx  <  m2  and  nx  <  n2.  The  total  cost  of  the  i-th  iteration  is  then  bounded  by 

C(in,  if)  +  c(2N)  log  (2V) 

<  C(2N ,  N )  +  c(2N)  log  (2N)  (3.56) 

<  C(2‘,  2i_1)  +  ci2‘ 

operations,  for  an  appropriate  constant  c.  Consequently,  we  have 

Theorem  3.6.  Given  that  0  ^  n  <  m  <  2n,  Algorithm  1  can  compute  the  scaled  Pade  fraction 
of  type  ( m,n )  for  A(z)  t  U  in  time  0 (n  log2  n). 


Proof:  Consider  the  recurrence  relation 


■ 
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C(  2k+\2k)  =  2[  C^,?"1)  +  ci  2‘],  k>  1, 

i- 1 

where  c  is  a  positive  constant.  Then 

C(2t+1,  2*)  =  X  [C(2^,  2'-1)  +  ci*]  +  C(2*,  2*"1)  +  ck2k 

i-  i 

=  2C(2*,  2*_1)  +  c*2*. 

With  n  =  2*,  results  on  recurrence  relations  (see  Bentley,  Hanken  and  Saxe  [££7V80])  then  yield 

C(2«,  n)  =  n[C(2,  1)  +  c£i] 

1-  1 

=  n[C(2,  1)  +  ck(k  +  1)  /  2  ] 

=  0(n  log2  n ). 

The  theorem  now  follows,  since  from  equation  (3.56)  for  m  and  n  satisfying  0  <  n  <  m  ^  2n 

C(m,n)<ii[C(2!,2l-1)  +  ci2!].  «* 

/-i 

Lemma  3.7:  Let  m  >  n,  and  let  A(z)  t  II.  Determine  an  integer  8  >  1  such  that 

A(z)=Al(z)  +  2--"  +  8A2(z),  (3.57) 

where 

At(z)  =  A(z)  mod  zm~n  +  1  (3.58) 

and  A2(0)  =£  0  if  8  <  oo.  If  8  <  n,  let  Snn-t(z)  /  Tn n_6(z)  be  the  scaled  Fade  fraction  of  type 
(n,  n-b)  for  1  M2(z)  .  Then  the  scaled  Pade  fraction  Smn(z)  /  Tmn(z )  of  type  ( m,n )  for  A(z)  is 
given  by 
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(•A1(z)S/1>fl-»(7)  +  if  8  <  n 


(3.59) 


otherwise, 


(S^-fiCz)  ,  if  8  <  n 


(3.60) 


otherwise. 


Proof:  If  8  <  n,  then 


a(S,„)  =  s(aa,.^  + 

^  max{  (m  -  n)  +  n,  (  m  —  n  +  8)  +  (n  -  8)} 


m . 


Also,  d(Tmn )  <  n.  Furthermore,  either  d(Smn )  =  m  and  d^T^)  =  n  or  both.  Moreover, 
A(z)  •  T^iz)  -  Smn(z) 


=  (Ax(z)  +  zm_'l+8A2(z)}  S„,*_8(z)  -  {A1(z)5n,n_6(z)  +  z'"-'’+er„,n_6(z)} 
=  zm_"+8  {A2(z)5„)fl.8(z)  -  r„,n_8(z)} 

=  zm-"+8{A2(z)0(z2"-8+1)} 

=  0(zm+n  +  l), 


Final]  y, 


GCD(Sm„Tm) 

=  GC£>(A,S....,  +  S,,.-.) 

=  GCD(z-— *7.,.-..  Si,.-.) 

=  GCD(Tn  n_8,  5n  n_8). 

If  8  >  n,  then  dearly  d(Sm„)  =£  m  and  3(7^)  =  n.  Moreover, 
A(z)  •  T^(z)  -  Sm„(z) 


=  zn{A1(z)  +  zm  n  +  8A2(z)}  —  znA1(z) 
=  zm+*A2(z) 

=  0(zm+n  +  1). 


L 
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Finally,  GCD(Smn,  T^)  =  z?  .  h 

Theorem  3.8:  For  arbitrary  m  >  n,  Algorithm  1  can  compute  the  scaled  Pade  fraction 
of  type  (m,  n )  for  A(z)  e  U  in  time  0(m  log  m)  +  0(n  log2  n). 

Proof:  If  8  >  n  in  equation  (3.57),  the  result  is  trivial. 

If  8  ^  n,  computation  of  the  first  2n  —  8  +  1  terms  of  1  /  A2(z)  requires  0(n  log  n) 
operations.  Using  (3.55), 

C{n,n~  8)  s  C(n,n), 

and  consequently  by  Theorem  3.6,  it  follows  that  the  cost  of  computing  the  scaled  Pade  fraction 
S*,*_8(z)  /  rn>„_8(z)  for  1/A2(z)  is  bounded  by  0(n  log2  n)  operations.  Finally,  the  cost  of 
computing  Smn(z)  in  equation  (3.59)  is  0(m  log  m)  operations.  a 


3.4.  Fast  Off-diagonal  Algorithm  for  a  Quotient  Povrer  Series 

Let  A  (z) ,  B  (z)  be  unit  power  series  and  let 

C(z)  =  -A(z)  /  B(z)  (3.59) 

be  the  quotient  power  series.  Given  the  non-negative  integers  m  >  n,  then  Algorithm  1  OFF- 
DIAG  can  be  used  to  compute  the  scaled  Pade  fractions  ^(z)  /  Tmn(z )  of  type  ( m,n )  for  C(z) .  As 
a  result, 

A(z)Tm,(z)  +  B(z)S„.(z)  =  0(z"*-1)  .  (3.60) 

Before  applying  the  algorithm,  the  quotient 

C(z)  mod  zm+"  +  1  =  -A{z)/  B{z)  mod  zm  +  n  +  1  (3.61) 


must  be  calculated.  However,  this  division  need  be  computed  modulo  zm+n+1,  only,  by  modifying 
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Algorithm  1  as  follows: 


ALGORITHM  2  :  OFFDIAG 
INPUT:  A(z),  B(z),  m,  n  ,  where 


(1)  m  and  n  are  non-negative  integers  with  m  >  n  ,  and 

(2)  A(z),  B{z)  are  unit  power  series  .  (  Note  that  only  A(z)  mod  zm+n  +  i  and  B(z)  mod  zm~',  +  1 
are  required ). 


OUTPUT: 


'  Si 
Tx 


S  o' 
T0 


,  where 


(1)  Sl  /  Tx  is  the  scaled  Pade  fraction  of  type  (m,  n)  for  -A(z)  /  fi(z),and 

(2)  S0  /  T0  is  the  predecessor  scaled  Pade  fraction  of  type  (m— \mn— 1,  for 

— A(z)  /  B(z),  given  that 

zm  =  GCD(S1,Tl)  . 


Step  1:  #  Initialization  # 


/  -  —  1 
M  -  (m  —  n) 
N  -  0 


S i  S0 

~A(z)/fi(z)  {mod  zM  +  1)  zM~1' 

n  r0 

♦- 

1  0 

Step  2:  #Cal culation  of  step-size  # 


i  -  i  +  1 

j  -  min  {2J  —  N  ,  n  —  N} 
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Step  3:  #  Termination  criterion  # 

If  5  =  0  then  exit 

Step  4:  #  Calculation  of  scaling  factor  for  Si  /  7\  # 

Determine  \MN  such  that 
zXmv=  GCD(  S\,Ti  ) 

Step  5:  #  Computation  of  residual  for  (z)  =  Sx  /  Tx  # 

Compute  and  7?!  such  that 

(AT,  +  BS,)  mod  =  z"""^****1 

where  tf,(0)  *  0  if  m,*  <  2 s  +  \m. 

Step  6:  #  Identification  of  Cases  # 

if  \^mn  —  s 

then  #  Case  of  Theorem  3.1  # 


Si  S0 

Si  So 

V  0 

Ti  r0 

Ti  T0 

0  1 

go  to  step  11 

else  #  Case  of  Theorem  3.4  and  3.5  # 
go  to  step  7 


#  Calculation  of  degrees  for  residual  scaled  Pade  fractions  # 


m  *-  s  +  \MN 
n  —  s  —  \iMN  —  1 


Step  7: 


■ 
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Step  8:  #  Computation  of  residual  for  =  S0  /  P0  # 


Compute  R0  such  that 


(A T0  +  B-S0)  mod  zM+N  +  m+*~7^  =  R0  i**"'2'*1'1 


where  tf0(0)  =£  0. 


Step  9:  #  Computation  of  residual  scaled  Pade  fractions  # 


Si  S0 
Pi  T0 


OFFDIAG  (  R0i  /?!,  m,  n  ) 


Step  10:  #  Advancement  of  scaled  Pade  fraction  computation  # 


Si 

So' 

*1 

Sol 

x  i 

To 

Pi 

i 

0 

2xmV-  i*  aw4' 


*i 

5o 

Pi 

Po 

. 

Step  11:  #  Calculation  of  degrees  of  Sx  /  Tx  # 


M  -  M  +  s 
N  -  N  +  s 
goto  step  2  ■ 


Algorithm  2  is  a  generalization  of  Algorithm  1,  since  it  can  be  used  to  produce  scaled  Pade 
fractions  for  a  single  power  series  A(z)  simply  be  setting  B{z)  =  —  1.  It  differs  from  Algorithm 
1  in  that  the  division 

A(z)  =  —  tf0(z)  / /?i(z)  mod  zm  +  "  +  1  (3.63) 

in  step  9  of  Algorithm  1  is  avoided.  Instead,  the  division  is  delayed  (i.e.,  immediately  subsequent 
to,  rather  than  prior  to,  the  recursive  call  of  OFFDIAG)  until  the  initialization 
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Sr  =  -A(z)/B(z)  mod  2m~n  +  l 


(3.64) 


in  step  1  of  Algorithm  2. 

There  are  also  various  sign  changes  introduced  in  steps  1,  5,  8  and  10.  These  account  for 


the  fact  that  the  algorithm  deals  with  the  power  series  ~A(z)  /  B(z)  rather  than  A(z)  /  B{z). 


This  notational  change  simplifies  the  development  of  subsequent  results. 

For  implementation  purposes,  Algorithm  2  can  result  in  considerable  savings  in  cost.  Practi¬ 
cally,  fast  division  is  significantly  slower  than  fast  multiplication,  by  an  asymptotic  constant  of 
approximately  7  (see  Verheijen  [VE7?83]).  For  example,  if  — A(z)//J(z)  is  normal,  then 
m  —  n  +  1  =  2  and  the  division  in  step  1  of  Algorithm  2  becomes  trivial.  However,  the  asymp¬ 
totic  cost  of  Algorithm  2  remains  the  same  as  the  asymptotic  cost  of  Algorithm  1  given  in  section 


3.3. 


The  proof  of  the  correctness  of  Algorithm  2  is  nearly  identical  to  the  proof  of  correctness  of 
Algorithm  1,  and  therefore  it  is  not  given. 

3.5.  Classical  Off-diagonal  Algorithm  for  a  Quotient  Power  Series 

Let  Si  /  Ti  be  the  scaled  Pade  fraction  of  type  (A///)  for  -A(z)/£(z)  such  that 


zXmn  =  GCD(Si  ,  TO  =  1. 


That  is,  \MN  =  0.  Then 


A  ■  Tx  +  B  •  Sx 


where  if  \iUN  <  <».  Thus,  (z*  Sx)  /  (z*  7\)  is  the  scaled  Pade  fraction  of  type 

(M+k  ,  N+k)  for  -A(z)/fl(z)  whenever  0  <  k  <  \iMN  .  If  \x.MN  <  »,  the  next  distinct  Pade  frac¬ 
tion  for  -A(z)  /  £(z)  along  the  (M-N)-th  off-diagonal  path  is  of  type  (M+s  ,  N+s),  where 


t 


' 
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S  =  P-m*  +  1  • 

By  so  selecting  the  step-size  s ,  Algorithm  2  may  be  used  to  compute  all  Pade  fractions  along 
the  ( M-N)-th  off-diagonal  path.  With  this  choice  of  s,  step  7  in  Algorithm  2  gives 

m  +  1 

n  =  0  . 

Consequently,  the  recursive  call  of  OFFDLAG  in  step  10  reduces  to 


Sv 

So 

—Rq/Ri  mod  z^  +  2 

Ti 

To 

1 

0 

and  step  11  then  yields 


'Si 

So 

Si-Si 

+  z 

•50 

z1*"'  ■  S, 

Ti 

To 

TvSi 

+  z"'^2 

To 

i™  ■  T ! 

Since  GCD{  Si  ,  7\  )  =  1,  it  follows  from  the  proof  of  theorem  3.4  that 

GCD  (  Si  •  Si  +  z*1^2  -S0  ,  Ti  •  Si  +  z^2  •  T0)  =  1  . 

Thus,  the  above  computations  can  be  repeated  for  the  scaled  Pade  fraction  of  type 
(M  +  s  ,  N  +  s). 

The  full  details  are  provided  in  Algorithm  3,  below.  To  simplify  the  presentation  of  subse¬ 
quent  results,  at  the  i-th  iteration,  the  scaled  Pade  fractions  Si  /  Tx  of  type  (M ,N)  is  denoted  by 
Si  /  Ti  and  its  predecessor  S0  /  T0  by  St-i  /  7i_i.  In  addition,  the  residual  power  series  and  R0  are 
denoted  by  Rt  and  Rt-i,  respectively. 


ALGORITHM  3:  OFFDIAG 


INPUT:  A,  B,  m,  n  ,  where 


a  ' 
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(1)  m  and  n  are  non-negative  integers  with  m  >  n  ,  and 

(2)  A  and  B  are  unit  power  series  .  (  Note  that  only  A  mod  zm  +  n  +  1  and  B  mod  zm+n+x  are 
required  ). 


OUTPUT: 


$i  + 1  Si 

Tl+ 1  T, 


where 


(1)  Si+X  /  Ti  +  1  is  the  scaled  Pade  fraction  of  type  (m,  n )  for  -A/B,  and 

(2)  St  /  Tj  is  the  scaled  Pade  fraction  of  type  n-X^-1)  for  -A/B,  given  that 

=  GCD(St+i,Ti  +  l)  . 


Step  1:  #  Initialization  # 


i - 1 

M  -  (m  —  n) 
N  -  0 


^+1 

5/ 

— A(z)/B(z)  mod  z^  +  1 

zM~1  ' 

r,+1 

71 

1 

0 

Step  2:  #  Termination  criterion  # 


If  N  =  n 
then  exit 
else  i  -  i  +  1 


Step  3:  #  Computation  of  residual  for  yMN  -  St  /  Tt  # 


Compute  \lmn  and  R,  such  that 


(A  Ti  +  B  Si )  mod  z 


M  +  N  +  2  3  _ 
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where  ^(0)  ^  0  if  ^iMN  <2  (n  -  N). 
Step  4:  #  Calculation  of  step- size  # 

5  -  min  {  \lmn  4-  1  ,  n-N  } 

Step  5:  #  Identification  of  Cases  # 

if  Mwv  —  s 

then  #  Case  of  Theorem  3.1  # 


Si  + 1 

Si  St-i 

V  0 

r1+i  r. 

Ti  T,.x 

0  1 

go  to  step  8 

else  #  Case  of  Theorem  3.4  and  3.5  # 
go  to  step  6 

Step  6:  #  Computation  of  residual  for  =  S/-i  /  7',- 1  # 

Compute  such  that 

(A  •  T/-X  +  B  ’  5^_i)  mod  z*  N  +  1  =  zM*N~1  , 
where  ^-1(0)  =£  0. 

Step  7:  #  Advancement  of  scaled  Pade  fraction  computation  # 


Si  + 1  Si 

Si 

Si- r 

1  0 

-Ri-.i  /  Ri  mod  ztl'*/+2 

2^ 

T,+i  ^ 

Tt 

r,-i 

0  z’**"*2 

1 

0 

Step  8:  #  Calculation  of  degrees  of  5i+i  /  7i+1  # 


N  ~  N  +  s 
M  —  M  +  s 


L 


■ 
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go  to  step  2  *> 


In  the  case  that  £(z)  — —  1  and  m  —  n,  Algorithm  3  is  precisely  the  algorithm  given  by  Cabay 
and  Kao  [CAS  83]  for  computing  diagonal  Pade  fractions  for  a  single  power  series.  Their  0(n 2) 
algorithm  (  since  steps  3  and  6  require  only  \xMN  +  2  terms  of  a  product  polynomials,  and  since 
one  of  the  multipliers  in  step  7  is  a  polynomial  of  degree  4-  1 ,  for  small  \iUN ,  no  advantage 
can  be  gained  from  using  fast  methods  for  polynomial  arithmetic)  is  shown  to  be  faster  than  other 
0(n 2)  algorithms  for  computing  diagonal  Pade  fractions,  such  as  those  of  Trench  [TRE6S]  and  Ris- 
sanen  [S/573]. 

The  more  interesting  observation  about  Algorithm  3  is  made  in  the  next  chapter,  however.  It 
is  shown  that,  when  the  given  quotient  power  series  is  a  rational  function,  Algorithm  3  along  one 
specific  off-diagonal  path  corresponds  exactly  to  Euclid’s  extended  algorithm  for  computing  the 
greatest  common  divisor  of  the  numerator  and  denominator  of  the  given  rational  function.  In  this 
sense,  Algorithm  3  is  a  generalization  of  Euclid’s  extended  algorithm.  It  is  for  this  reason  that  we 
choose  to  call  Algorithm  3  classical . 


CHAPTER  4 


GREATEST  COMMON  DIVISOR  COMPUTATIONS  OF  POLYNOMIALS 

4.1.  Introduction 

In  this  chapter,  Algorithm  2  and  Algorithm  3  are  examined  when  they  are  applied  to 
— A(z)/B(z),  where  A(z)  and  £(z)  are  polynomials  of  degrees  m  and  n  respectively.  In  addition,  it 
is  assumed  that  A(0)^0  and  5(0)^0. 

To  permit  the  analysis,  we  need  the  following 
Definition:  The  reciprocal  of  a  polynomial 


P(z)  =  Po  +  Pi  Z  +  '  ■  •  +  Pn? 


(4.1) 


of  at  most  degree  n  is  defined  to  be 


p*(z)  =  p0  z"  +  Pi  z"  1  +  •  •  •  +  pn  =  z*  P{z  :). 


(4.2) 


The  name  originates  from  the  fact  that  the  zeros  of  P*(z)  are  the  reciprocals  of  those  of  P{z)  if 

PoPn  *  0. 


Clearly 


[p*  wr  =  P(?)  ; 


(4.3) 


that  is,  the  operation  of  forming  the  reciprocal  polynomial  is  involuntary.  Also 


[P(z)  •  Q(z)Y  =  P«  (z)  •  Q*  (z), 


(4.4) 


and 


[z*  P(z)Y  =  P*  (z). 


(4.5) 
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Let  A(z)  be  a  unit  power  series.  The  first  n  +  1  terms  of  A(z)  are  denoted  by 

4.  00  =  2  ^  ^ 

<-o 

=  A(z)  mod  z"+1. 

The  reciprocal  A*(z)  of  A*  (z)  is  then  given  by 

A?(z)  =  ^  ■ 
i-0 


(4.6) 


(4.7) 


4.2.  Eudid’s  Extended  Algorithm 

Theorem  4.1.  Algorithm  3  applied  to  — A(z)/B(z),  where  A(z)  and  B(z )  are  polynomials  of 
degree  m  and  n  respectively,  is  equivalent  to  Euclid’s  extended  algorithm  for  computing  the 
greatest  common  divisior  of  A*(z)  and  £*(z). 

Proof:  For  completeness  in  presentation,  assume  that  modulo  operations  are  not  performed  in 
steps  3  and  6  of  Algorithm  3.  Thus, 

AT,  +  B  S,  =  z"*"*  “"'*1  R,  (4.8) 

+  (4.9) 

where  R,(z)  and  /?/-i(z)  are  polynomials  of  degrees  n-N-  -1  and  n-N,  respectively  1  .  By 
taking  reciprocals  of  (4.8)  and  (4.9),  it  follows  that 

A*  •  7f  +  B?  •  Sf  =  /?f  (4.10) 

and 

AR  •  7*_!+  BR  ■  Sf_!  =  /#_!.  (4.11) 


:By  convention,  a  polynomial  of  negative  degree  is  the  zero  polynomial. 
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Let 


Qi+i(z)  =  /  Rfa)  fn°d  z*MS  +2 

be  a  polynomial  of  degree  jv  + 1  ■  Then 

fl/-i(z)  “  <2/+i(z)  •  Rt(z)  =  Z*MN  2  R(z), 
for  some  polynomial  /?(z)  of  degree  m—N—  \iMN  -2,  Taking  reciprocals  of  (4.13), 

r,'-iW  -  aii  w  •  woo  =  **w. 

That  is,  Q’+i(z)  and  7r(z)  are  the  quotient  and  remainder,  respectively,  c 

W-i(0  by  Rf(z). 

With  0i+1(z)  defined  by  (4.12),  step  7  of  Algorithm  3  yields 


Si  +  i  Si 

_|4A/N  +2  .  r  _  /O 

z  *J/-1  VZi  +  l 

$ 

z*MN  •  5/ 

Ti+ 1  7) 

_H-A/7V  +2  .  r  —  O 

z  ■*  i  - 1  tfi+1 

Ti 

which  on  taking  reciprocals  becomes 


i'l  if 

re*  _  /-t*  .  cR  c*i 

*>/-l  t/i+i 

Tf+i  7f 

it.  -  ef.i  •  7?  7? 

Looking  ahead  one  iteration,  we  obtain 

Kf,i  =  A'  ■  7*  ,  +  W  •  S’rl 

=  A"  ■  [7?.,  -  e?rt  •  7f]  +B"  ■  [  Sf_i  -  OS-t  '  5,'] 

=*f-i  -  a**i  • W- 

Thus,  W(z)  =  7?f+i(0  in  equation  (4.14). 

Summarizing,  let  Qf.i  be  the  quotient  on  division  of  l?f- ,  by  Rf.  Then 

nR  —  nR  _  rJt  .  nR 

*d  +  l  —  Ki-\  VzJ  +  1 

sf.i  =  s’. ,  -  e*  i  ■  if 


(4.12) 


(4.13) 
it  follows  that 

(4.14) 
n  division  of 


(4.15) 


(4.16) 


(4.17) 


(4.13) 
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tf+i  =  Tf-i  -Qf+i-T?, 

which  are  the  fundamental  relations  describing  Euclid’s  extended  algorithm. 

To  complete  the  proof,  it  remains  to  show  that  the  initial  conditions  for  Euclid’s  extended 
algorithm  are  satisfied.  Initialization  in  step  1  of  Algorithm  3  yields 

(4.19) 

where,  as  in  (4.12),  (2o  is  the  quotient  on  division  of  A*  by  BR .  Steps  3  and  5  with  i  =  0  then 
become 


50  5-1 

-Go 

1 

i _ 

1  0 

A  ■  1  +B  •  (  -Q0)  = 


(4.20) 


and 


A  ■  0  +  B  ■  =  r»-«- 1  R^ , 


(4.21) 


Taking  reciproicals  of  (4.19),  (4,20)  and  (4.21),  it  follows  that 


R*Li  =  B*  S*i  =  1  7^i  =  0 

Ro  =  A*  -  Q*  •  B*  S*0  =  -Q5  T*0  =  1  • 


(4.22) 


Corollary  4.2.  At  the  i-th  iteration  of  Algorithm  3,  let  St  /  T,  be  the  scaled  Pade  fraction  of  type 
(MAO  f°r  -A  /  B,  such  that 


A  •  T,  +  B  -  5/  = 


ZM+N+  +  1  £ 


Then, 


R,-l'S,  -  z’1"'*2  R,  ■  S,-1  =  (-1)'  -  A 

R, -i  ■  T,  -  z""*2  R,  ■  T,-i  =  (-1)'*1  •  B  (4.23) 

S, -i  T,  -  S,  ■  T,-!  =  (-I)'*1  z"*»-‘  . 


Furthermore,  the  algorithm  terminates  for  some  i  =  k,  where  k<n.  On  termination,  Sk+l  /  Tk+l  is 


I 
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the  scaled  Pade  fraction  of  type  ( m,n )  for  —A  /  B  such  that 


A  •  Tk  +  l  +  B  ■  Sk+i  —  0. 


In  addition, 


A  ■  Tk  +  B  ■  Sk 


m+n-2  A_  - 1 
Z 


Rk 

*  t 


where 


=  GCD(st ,  r,.,) 

and 


(4.24) 


(4.25) 


/?*  =  GCD(A  ,  B). 

Proof:  From  Euclid’s  extended  algorithm  (4.18)  with  initial  conditions  (4.22),  it  follows  that 

(see,  for  example,  McEliece  [Mc£78]) 

tff-1  -S*  -  Rf  '  Sf_l  =  (-1  )'  -  A* 

Rf-i  -7?  -  Rf  ■  7f_!  =  (-l)i+1  •  B*  (4.26) 

S*_i  •  7f  -  *7f.x  =  (-1)/+1. 

By  using  the  correspondence  established  in  theorem  4.1  and  taking  reciprocals,  equations  (4.2 6) 
result  in  equations  (4.23). 

Furthermore,  Euclid’s  extended  algorithm  terminates  for  some  k  <  n  when 

A*  •  7?+1  +  B*  •  Sf+1  =R?,i=0  (4.27) 

and 

A*  •  Tf  +  B*  ■  Sk  =  Rk  =  GCD  (A*  ,  B*),  (4.28) 

where  X  —  d(Bf),  3(5*  +  !)  ^  m  —  X,  d(7f+1)  <  n  —  X,  d(S£)  <  m  —  X  and  d(7f)  <  n  -  X. 


Taking  reciprocals  of  (4.28)  with  respect  to  zm+"  x  1  results  in 
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A  ■  Tk  +  B  •  Sk  =  zT+n~  x  _1  •  Rk, 

where  i?t(0)^0.  Thus,  Sk  /  Tk  is  the  scaled  Pade  fraction  of  type  (m-  X  -1  ,  n-  X  -1)  for 
-A  /  B.  Taking  reciprocals  of  (4.27)  with  respect  to  zm+"  gives 

A  ’  Tk+ 1  +  B  •  Sk+1  =  0,  (4.29) 

where  d(St+1)  ^  m  and  <9(ri+1)  —  n.  Thus,  Sk+1  /  Tk+l  is  the  scaled  Pade  fraction  of  type  (m,n) 
such  that 

GCD(Sk  +  1  ,  Tk+1)  =  /  .  (4.30) 

Thus  X=Xmn.  ■ 

As  a  consequence  of  Corollary  4.2.  Algorithm  3  computes  co-multipliers  Sk  and  Tk,  only, 
such  that 

A*  •  7?  +  B*  •  S'?  =  RRk  =  GCD(Ar  ,  BR)  .  (4.31) 

The  remainder  Rk  is  available  only  if  the  multiplications  in  steps  3  and  6  are  performed  without  the 
modulo  operation. 


4.3.  Fast  GCD  Computations 

Theorem  4.3.  Algorithm  2  applied  to  -A(z)  /  B{z),  where  A(z)  and  B(z)  are  polynomials 
of  degree  m  and  n,  respectively,  returns  the  scaled  Pade  fraction  Si  /  7\  of  type  ( m,n )  such  that 

A*  •  7?  +  BR  •  SRi  =  0  (4.32) 

and  its  predecessor  S0  /  T0  of  type  (m-  Xm  -1  ,  n  —  Xm„  -1)  such  that 

Ar  ■  T%  +  BR  ■  S$  =  ,  (4.33) 


where 
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**“  =  CCD(St  ,  T,) 
and 


P0  =  GCD(A  ,  B). 

Proof:  Since  scaled  Pade  fractions  are  unique,  the  result  is  an  immediate  consequence  of 

Corollary  4.2.  ■ 

The  greatest  common  divisor  R0  is  not  explicitly  computed  by  Algorithm  2.  However, 
d(R0)  =  \mn  <  n.  Using  fast  multiplication,  it  can  therefore  be  determined  in  0(n  log  n).  As  a 
consequence  of  this  and  Theorem  3.8,  we  have 

Theorem  4.4.  Algorithm  2  can  compute  the  greatest  common  divisor,  the  cofactors  and  the 
comultipliers  of  two  polynomials  of  degrees  m  and  n,  where  m  >  n,  in  time 
0(m  log  m )  +  0{n  log2  n).  ■ 

Thus,  Algorithm  2  for  GCD  computations  is  basically  of  the  same  asymptotic  complexity  as 
the  fast  algorithm  of  Moenck  [MOE13],  Aho,  Hopcroft  and  Ullman  [AH 07 4]  and  Brent,  Gustavson 
and  Yun  [PPfSO],  which  are  of  complexity  0(  (m+n)  log2  (m+n)  ).  However,  Algorithm  2  has 
the  advantage  of  being  partly  iterative  (approximately  half  as  many  recursive  calls  are  used  in  com¬ 
parison  with  the  other  fast  methods).  This  can  result  in  significant  cost  savings  in  an  implementa¬ 
tion  environment  (Verheijen  [V£P83]). 


4.4.  Antidiagonal  Computations 
Let 


A(z )  = 


2 

J-0 


(4.34) 


be  a  unit  power  series,  and  let 


57 


A*d(z)  =  £  ad_ t  zf  ,  0,  (4.35) 

;-o 

be  the  reciprocal  of  the  first  d+1  terras  of  A(z).  In  this  section,  we  examine  the  intermediate 
results  obtained  by  Algorithm  2  2  while  computing  the  scaled  Pade  fraction  of  type  (n,n),  2n<d, 
for  the  two  polynomials  1  +  z  A*  (z)  and  B*(z)  =  -1. 

At  the  i-th  iteration  of  Algorithm  2,  the  scaled  Pade  fraction  SNN{z)  /  TNN(z)  of  type  (N.N) 
for  -{1+z  A*(z)}  /  Z?*(z)  and  its  predecessor  SN*N*(z)  /  TN»N*(z)  are  determined  such  that 

{  1+z  Aj(z)}  •  Tm(z)  -  Sm(z)  =  z“ ' 1  R:(z)  (4.36) 

(1+  z  A?(z)}  •  T„v.(z)  -  S„v(z)  =  z2"'  +1  R„(z)  ,  (4.37) 

where 

W  =W-Xra-1 
and 

z'*"  =  GCD(SmTm)  . 

By  taking  reciprocals  of  (4.36)  and  (4.37),  it  follows  that 

A„(z)  •  7£v(z)  -  tf?(z)  =  z^1  {5jS„(z)  -  (4.38) 

(z)  •  (z)  —  ^?o(z)  =  z^  1  {5$^.  (z)  —  i*»N»(z)}  (4.39) 

where 

d(R! )  =  d  —  N  —  ilnn  (4.40) 

and 

d(/?S)  =  d  -  .  (4.41) 


2  It  is  assumed  that  no  truncation  of  polynomial  operations  takes  place. 


JVt  ‘ 


' 
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Theorem  4.5:  Let  M—d—N  and  M*  =  d  -  NT .  Then  the  scaled  Pade  fractions  of  type  (M  ,  N ) 
and  (M*  ,  N*)  for  A(z)  are  /?f  (z)  /  T%N(z)  and  R*(z)  /  ^.^.(z)  ,  respectively. 


Proof:  From  (4.38),  (4.39),  (4.40)  and  (4.41),  clearly  the  degree  and  order  conditions  for  scaled 
Pade  fractions  are  satisfied.  To  complete  the  proof,  let 

G(z)  =  GCD(7^  ,*?). 

Then,  from  (4.38),  G(z)  must  divide  not  only  T$N  but  zd+l  S*N(z),  as  well.  Since 
GCD(S*n  ,  7^*)  =  1,  it  follows  that 

G(z)  =  zx 

for  some  \  ^  0.  Thus,  R* (z)  /  T%jy(z)  is  a  scaled  Pade  fraction  of  type  (M^V).  Similarly, 
/?o(z)  /  Tn*n*(z)  is  a  scaled  Pade  fraction  of  type  (A/* ,  N*).  ■ 

Observe  that 


M  +  N  =  d 


(4.42) 


and 

M*  +  N*  =  d  .  (4.43) 

The  scaled  Pade  fractions  /?? (z)  /  T%N(z)  and  /?o(z)  /  7^*  *(z)  both  lie  along  the  J-th  anti-diagonal 
path  of  the  scaled  Pade  table  T(A),  where 


Definition  (see,  for  example,  Brent  [B^ZtSO]).  The  <f-th  anti-diagonal  path  of  the  scaled  Pade 
table  T(A)  for  a  unit  power  series  A{z)  is  defined  to  the  set  of  all  scaled  Pade  fractions  of  type 
(i  J),  where 

i  +  j  =  d  .  (4.44) 

But,  by  Theorem  4.5,  the  scaled  Pade  fraction  (z)  /  T%N(z)  and  its  predecessor 
tf$(z)/^v00,  along  the  <f-th  anti-diagonal  path,  are  obtained  by  applying  Euclid’s  extended 
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algorithm  to  the  reciprocals  of  1  +  z  Ad(z)  and  BR(z )  =  —  1,  that  is,  to  the  polynomials 
Ad(z)  +  zd  +  l  and  —zd+1.  However, 

GCD(Ad(z )  +  /  +  1  ,  -zd+1  )  =  GCD{  Ad(z)  ,  -z^1  ).  (4.45) 

This  is  precisely  the  result  given  by  McElice  [McElS]]  namely,  Euclid’s  extended  algorithm  applied 
to  Ad(z)  and  -zd+l  yields  all  the  Pade  fractions  along  the  J-th  anti-diagonal  path  of  the  Pade  table 
for  A(z).  Equivalently,  by  Theorem  4.5,  the  same  Pade  fractions  can  be  obtained  by  applying 
Algorithm  3  to  the  two  polynomials  1  +  z  Ad(z)  and  BR{z)  =  —  1,  and  then  taking  reciprocals  of 
the  results. 

In  addition,  by  Theorem  4.5,  Algorithm  2  can  be  used  to  compute  any  specific  scaled  Pade 
fraction  R*  (z)  /  TRN(z)  of  type  (M  ,  N),  where  M  +  N  -  d,  along  the  d-\h  anti-diagonal  path  for 
A(z).  This  requires  0(  N  log2  N)  arithmetic  operations  to  compute,  by  Algorithm  2,  the  scaled 
Pade  fraction  5,v,v  (z )  /  Tm(z)  of  type  (iV,A0  for  1  +  z  Ad(z),  plus  an  additional 
0(  ( d—N )  log  ( d-N )  )  operations  to  determine  wiiich  satisfies  (4.36).  This  asymptotic  cost  is 
the  same  as  the  cost  of  the  fast  algorithm  of  Brent  [5^80]  for  computing  the  greatest  common 
divisor  of  Ad{z )  and  —  zd+1.  However,  their  PRSDC  algorithm  can  compute  a  specific  Pade  frac¬ 
tion  along  the  d-th  anti-diagonal  only  with  significantly  more  complications. 


CHAPTER  5 


THE  INVERSE  OF  HANKEL  AND  TCEPLITZ  MATRICES 


5.1.  Introduction 


In  this  chapter,  scaled  Pade  fractions  are  used  to  construct  the  inverse  of  a  nonsingular 
Hank  el  matrix 


Hnn  = 


a\ 


a „ 


a2n-\ 


and  the  inverse  of  the  Toeplitz  matrix 


aln-\ 


associated  with  the  Hankel  matrices. 


The  formula  derived  for  the  inverse  is  similar  to  those  given  by  Trench  [7RE65],  Zohar 
[Z0//74],  and  Kailath  et  al  [AA/78],  but  it  has  the  advantage  of  being  successful  in  the  degenerate 
case  when  a  principal  minor  of  Hnn  is  zero.  The  formula  given  by  Gohberg  and  Semencul  [£/?£80] 
does  succeed  for  the  degenerate  case,  however,  it  is  a  slower,  0(n2),  algorithm. 


5.2.  The  Inverse  Formula 
Let 
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S(z)/T{z)=  2>zV  2^  (5.1) 

1-0  1-0 

be  the  scaled  Pade  fraction  of  type  (n,n)  for 

oc 

Mz)=2qz'-  (5.2) 

*-0 

Assume  that 

a*  =  0  ■  (5.3) 

Since  $(z)/T(z)  is  also  a  Pade  form  of  type  (n,n),  the  denominator  T(z)  satisfies 


ax 

&n  +1 

1 - 

vT  • 

1 

'  0  ' 

an- 1  . 

a* 

a2n-2  Clin- 1 

a2n-l  0 

h 

h 

— 

.  o 

_ i 

Furthermore,  from  Collorary  2.8,  det(//M)  ^  0  if  only  if  t0  0.  In  the  remainder  of  this 
chapter,  assume  det(H^)  ^  0.  We  now  shall  provide  an  algorithm  for  computing  the  inverse  of 

hm. 

Let 

U(z)  /  V(z)  =  £  ut  <  2  viz>  (5.5) 

/-0  i-0 

be  the  scaled  Pade  fraction  of  type  (n-l,n-l)  for  A(z)  (the  predecessor  of  S(z)  /  T(z)  ).  Then 
V(z)  satisfies 
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v«- 1 

1 

al 

an 

• 

0 

• 

* 

• 

•  •  •  a2n-2 

Vl 

0 

V0 

Lemma  5.1.  Let 


•y  =  On>  •  •  •  ,a2i.-l)-(v»-l,  •  ■  •  ,V0)‘. 


Then  y  ^  0. 


(5.6) 


(5.7) 


Proof:  Suppose  7  =  0.  Then,  from  (5.6),  (vn_l5  •  •  •  ,v0,0)  *  satisfies 


a: 

an  tf«  +  l 

v*-i 

■  0 

<*n-l  ■  ■  • 

*n 

a2i\-2  a2n-\ 

a2n-l  0 

v0 

0 

= 

0 

that  is,  (vn_i  ,  •  •  •  ,v0  ,  0)'  is  also  a  solution  of  (5.4).  This  implies  det {HM)  =  0,  which  is  a  con¬ 
tradiction.  ■ 


Lemma  5.2:  Let  —  (^)u-i.  denote  the  inverse  of  HM.  Then 

blA  =  Kti  =  v*_,  /  7  ,i=  1,  ■  ■  ■  ,  n, 

where  7  is  given  by  (5.7). 


Together  with  7,  system  (5.6)  becomes 


a\ 

an 

v„-i 

•  0 

'  '  i* 

•  ■  ■  a2n-2 

•  •  •  a2n-l 

v0 

II 

0 

.  7 

(5.8) 


Proof: 


r 


63 


Equivalently, 


*  1 

_ 1 

- 

Vo 

'u 


bln  0 


Kl  ■ 


0 

y 


and  (5.8)  now  follows  since  BM  is  symmetric.  h 

To  obtain  the  remaining  elements  of  Bnn,  two  cases,  tn  =£  0  and  r„  =  0,  of  the  solution  of 
(5.4)  are  identified. 


Case  1:  Assume  that  tn  #  0.  Then  (5.4)  becomes 


«i 

• 

+ 

an-\ 

an 

a2 


On 

+ 1 


°n  an  + 1 


°2a-2  a^-i 
a2n-l  0 


■ 

1 

o 

» 

h 

— 

•  o 

h 

o 

(5.9) 


Let 


*2 


a. 


-*«  +  ! 


an  a, 


•d+i 


a2n-2  0-2n-\ 

a2*-\  0 


(5.10) 


Then  det(//n  +  1  „)  =£  0  ,  for  otherwise  (5.9)  also  has  a  solution  with  tn  =  0.  This  is  a  contradiction 
of  the  uniqueness  (up  to  a  scalar  )  of  the  solution  of  (5.4). 

Lemma  5.3:  Assume  tn  ±  0,  and  let  Bn+l  n  =  (fy/)Jj-i  denote  the  inverse  of  Hn  +  l  n. 


Then 


.*  r 


. i 


ij  =  1,  '  •  *  , 


(5.11) 


bij  —  b,^ij  (  +  /  to)  bnj  , 

where  b0j  =  0,  j  =  1,  •  •  •  ,  n. 


Proof:  From  (5.4) 
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b ii  ‘  •  bu 

an  +  1 

- 

[K 1,  •  ■  •  ,bM]  =  ~t0 

. 

; 

• 
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aln-l 

h 

.  *«1  ‘  ’  bM  _ 

0 

\bn  i >  '  ‘  '  j  bn 


However, 


O 

o 

1 

an  +  l 

*n  bu 

• 

Hnn 

■ 

=  In~ 

a2n-\ 

bn- 1,1 

0 

Consequently,  by  substituting  (5.13)  into  (5.12),  it  follows  that 


’  *11 

•  .  . 

bm  ' 

- 1 

O 
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• 

• 

_ 

*11 

...  *1(1 

bnl 

.  .  . 

bnn 

•  •  •  bn  -l,n 

■to1 


tn  bn\  ‘  tn  bM 


t\  bn  i  •  •  •  fj  bn 


]  .  (5.12) 


(5.13) 


(5.14) 


and  (5.11)  follows.  ■ 
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Case  2.  Assume  tn  =  0.  Consider  the  subsystem 


a2 


On 


an 


a2*-2 


f*-l 

• 

—  ~*0 

• 

h 

a2/i-l 

(5.15) 


of  (5.4).  Then  the  matrix  on  the  left-hand-side  of  (5.15)  is  nonsingular.  For,  suppose  otherwise. 
Then  there  exists  a  nontrivial  vector  [  an_1?  •  ■  ■  ,0^]'  such  that 


a2  • 

an 

a«-l 

o 

‘ 

• 

• 

’ 

= 

• 

a*  • 

a2n-2 

•  a 

0 

However,  by  taking  the  transpose  of  (5.15),  it  follows  from  (5.16)  that 


a2  •  •  *  an 

1 

[  £*  + 1)  '  '  '  >a2*-l] 

• 

—  1  ‘  )  C] 

• 

8  • 

an  '  <*2a-2 

«1 
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0 


(5.16) 


(5.17) 


Consequently, 


a\ 


On 


an 


a2n-\ 


5 


which  implies  that  det  {HM)  ±  0.  This  contradiction  implies  that 


. 


. 
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det 


a2  •  *  '  a. 


an  '  •  •  a2jt_ 2 


*  0. 


Now,  since  tn  =  0,  it  follows  from  (5.4)  that  det  (//n+1  „)  =  0.  Let 


= 


«2 


+ 1 


<3- 


fln  +  l 


a2n-2  £*2*- 1 

a2*-l  1 


(5.18) 


Then,  expansion  with  respect  to  the  last  column  yields 


det(//*+ln)  =  det 


a2  ■  ■  ■  an 


an  +  1 


a2n-2 


*  0. 


(5.19) 


Lemma  5.4.  Assume  tn  =  0,  and  let  =  {b*j)h- i  <^enote  the  inverse  of  /£+ltlI.  Then 

bij  =  +  (bi*-tn  +  i-t  /  to)  h*j  ,/,  7  =  1,  ■  ■  '  ,n,  (5.20) 

where  b*0j  =  0,  7  —  1,  •  •  •  ,  n  . 

Proof:  The  result  follows  using  arguments  similar  to  those  of  the  proof  of  Lemma  5.3.  ■ 

Theorem  5.5  Let  (f„,f*_  1,  ■  •  •  ,t0)‘  and  (vrt_x,v„_2,  •  •  ■  ,v0)‘  be  given  by  (5.1)  and  (5.5), 


respectively.  If  r0  =£  0,  then  the  inverse  of  the  HM  is  given  by  the  Christoff el-Dar bo ux  formula 


I  ■ 
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tn 


(5.21) 


Proof:  Consider  the  case  tn  =£  0.  Then  equation  (5.11)  of  Lemma  5.3  yields 

bM  =  ~  (t0  /  t„)  bu 


and 


bi+i ,n  =  bi  n  -  (r„_,  /  r0)  K 


=  1,  2, 


n , 


Consequently, 

Vi,*  =  bi<n  +  (f„_,  /  tn)  bu  ,  t  =  1,  2,  ■  ■  ■  ,n.  (5.22) 

Since  ^  ti  ,  equation  (5.22)  inserted  into  (5.11)  gives 

by  =  b^ij  —  /  fo)  [  bj+ 1,„  —  (r„_;  /  f*)  ]  . 

Thus, 

bj+i,i-\  —  bji-i  —  (tn-j  /  t0)  [  bj  „  —  (t„  +  1_i  /  f„)  i>i„  ]. 

Therefore, 

bij  =  bj+i'i-i  +  t0  1  [  tn-j  bi  n  —  tn  +  i_i  bj+i'H  ]. 

From  (5.8)  of  Lemma  5.2,  it  now  follows  that 

bij  =  bJ+u  + 1  +  (y  f0)-1  [  tn-j  v„_!  -  tn+1-t  v„_w  ], 
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i,  j  =  1,  2,  •  •  •  ,  n,  (5.23) 

where  all  variables  with  subscripts  outside  the  range  1,  2,  •  •  •  ,  n  are  assumed  to  be  zero.  But 
equation  (5.23)  corresponds  exactly  to  (5.21). 

The  case  tn  =  0  makes  use  of  Lemma  5.4,  and  the  proof  proceeds  in  a  similar  manner.  ■ 


Now  consider  the  nonsingular  Hankel  system 


Hn<nx  = 


a2n-l 


x  —  b 


(5.24) 


for  an  arbitrary  vector  b.  From  (5.21),  the  computation  of  x  =  H •  b  can  be  done  by  per¬ 
forming  four  polynomial  multiplications.  Using  fast  polynomial  multiplication  methods,  this 
requires  0(n  logn  )  arithmetic  operations.  By  theorem  3.6,  the  polynomials  T(z)  and  V(z )  can  be 
computed  in  0(n  log2  n )  operations,  and  therefore 


Theorem  5.6.  The  system  (5.24)  can  be  solved  m  0(n  log2  n)  arithmetic  operations. 


Example:  Let 


2  0  0 
0  0  1 
0  1  4 


J 


and  let 


A(z)  -flo  +  2z  +  z4+4z3  +  0/  + 

be  the  associated  unit  power  series.  Then  the  scaled  Pade  fraction  of  type  (3,3)  for  A(z)  is  given 
by  5(z)  /  T(z),  where 

S(z)  =  2  a0  +  4  (1  —  2 a0)  -  16  (1  -  2a0)  z2  +  (64  -  a0)  z3 


- 


' 
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and 

7(z)  =  2  -  8z  +  32Z2  -  z3  . 

Thus,  //3  3  is  nonsingular,  since  7(0)  =  t0  #  0.  The  predecessor  of  S(z)  /  7(z)  is  the  scaled  Pade 
fraction  t/(z)  /  V(z)  of  type  (2,2),  where 

U(z)  =  a0  z  +  2  z2 

and 


V(z)  =  z. 


Therefore,  the  solution  of  system  (5.4), 
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is  (^^2dido)‘  =  (-1,  32,  -8,  2)';  and  the  solution  of  system  (5.6), 
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O' 
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Vl 

0 

Vo 


is  (v2,  vl5  v0)‘  =  (0,  1,  0). 


The  inverse  H3  3  is  obtained  by  the  Christoffel-  Darboux  formula  (5.21)  to  be 
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Since  Toeplitz  matrices  are  equivalent  to  Hankel  matrices  (simply  by  a  change  of  notation), 
the  results  of  theorem  5.6  apply  to  Toeplitz  matrices,  as  well. 

Consequently,  by  Theorem  5.6,  the  inverse  of  a  Toeplitz  matrix  of  order  n  can  be  obtained  in 
0(n  log2  n)  arithmetic  operations.  This  is  done  with  one  call  of  Algorithm  1  for  the  computation 
of  a  specific  scaled  Pade  fraction.  The  fast  algorithm  of  Brent  [5&E80]  is  also  of  complexity 
0(n  log2  n),  but  the  computation  of  the  Toeplitz  inverse  may  require  as  many  as  four  calls  of  their 
0{n  log2  n )  PRSDC  algorithm  for  computing  a  specific  Pade  fraction.  In  addition,  each  call  is  bur¬ 
dened  with  complications,  alluded  to  earlier  in  Section  4.4. 


CHAPTER  6 


CONCLUSIONS  AND  SUGGESTIONS  FOR  FUTURE  RESEARCH 


Central  to  the  classical  theory  of  Pade  approximants  of  power  series  are  the  concepts  of  the 
Pade  form,  which  always  exists  but  may  not  be  unique,  and  the  Fade  fraction,  which  is  unique  but 
in  a  certain  sense  may  not  exist.  The  fundamental  definition  introduced  in  this  thesis  is  the  scaled 
Pade  fraction  which  exists  uniquely.  It  is  shown  that  scaled  Pade  fractions  satisfy  a  three-term  rela¬ 
tionship  between  elements  along  an  off-diagonal  path  of  the  scaled  Pade  table.  This  relationship 
circumvents  the  problems  of  the  degenerate  case  -  i.e.,  the  problems  which  plague  other  relation¬ 
ships  such  as  those  upon  which  the  e-algorithm  and  the  qd-algorithm  are  based. 

The  three-term  relationship  is  used  to  develop  an  (^(n2)  algorithm  (  Algorithm  3  )  which 
computes  along  an  off-diagonal  path,  a  finite  sequence  of  successive  scaled  Pade  fractions  for  the 
quotient  of  two  power  series.  In  the  case  that  the  power  series  is  normal,  Algorithm  3  is  identical 
to  the  one  described  by  Brezinski  [BRZ16].  In  the  case  that  computations  progress  along  the  diago¬ 
nal,  Algorithm  3  becomes  that  of  Cabay  and  Kao  [CAB  83].  Furthermore,  it  is  shown  that  if  the 
two  power  series  are  finite  (i.e.,  polynomials),  then  Algorithm  3  with  computations  along  one 
specific  off-diagonal  path  is  exactly  equivalent  to  Euclid’s  extended  algorithm  for  computing  the 
greatest  common  divisor  of  two  polynomials.  However,  other  off-diagonal  paths  can  be  used  to 
compute  greatest  common  divisors,  and  it  remains  a  subject  for  future  research  to  determine  the 
optimal  one. 

By  doubling  the  step  at  each  iteration,  the  three-term  relationship  gives  rise  to  an 
0(n  log2  n)  algorithm  (Algorithm  1  or  Algorithm  2).  The  cost  complexity  assumes  the  existence 
of  fast  methods  for  polynomial  arithmetic.  The  decision  to  double  the  step-size  (rather  than  to  tri¬ 
ple  it,  for  example)  is  not  accidental.  We  believe  that  this  choice  is  optimal  (within  the  constraints 
placed  by  fast  polynomial  arithmetic  methods),  but  a  formal  proof  remains  to  be  obtained. 
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When  applied  to  polynomials,  Algoritlim  2  at  each  iteration  routinely  produces  intermediate 
polynomial  remainder  sequences.  For  this  reason,  and  because  Algorithm  2  is  simply  a  fast  version 
of  Algoritlim  3,  Algorithm  2  is  truly  a  fast  Euclid’s  extended  algorithm.  The  PGCD,  HGCD  and 
EMGCD  algorithms  of  Moenck  [A/CZT73],  cf  Alio,  Hopcroft  and  Ullman  [AH 01  A],  and  of  Brent, 
Gustavs  on  and  Yun  [5/IE80],  respectively,  also  compute  greatest  common  divisors  with  the  same 
cost  complexity,  but  these  methods  are  simply  GCD  methods.  They  do  not  produce  intermediate 
polynomial  remainder  sequences.  In  addition,  the  artificial  splitting  cf  polynomials  used  by  these 
methods  to  achieve  their  speed  makes  them  difficult  to  understand.  Finally,  these  methods  are 
totally  recursive,  which  makes  them  more  costly  to  implement  than  Algorithm  2  which  is  semi- 
iterative. 

For  computing  scaled  Pade  fractions  for  a  power  series,  Brent  et  al  [iWJZsSO]  have  shewn  that 
the  EMGCD  algorithm  can  be  modified,  with  substantial  extra  detail,  to  compute  entries  along  an 
anti-diagonal  path  of  the  Pade  table.  We  have  shown  that  Algorithm  2  and  Algorithm  3  can  also 
be  used  to  compute  such  entries,  routinely.  From  a  practical  point  of  view,  however,  it  seems  to 
us  that  computations  along  an  anti-diagonal  path  are  not  as  natural  as  they  are  along  an  off- 
diagonal  path.  If  the  choice  of  the  anti-diagonal  path  is  incorrect,  computations  must  be  restarted. 
For  off-diagonal  computations,  n  need  not  be  known  apriori,  and  from  an  application  point  of 
view,  this  may  be  one  of  the  most  important  contributions  of  this  thesis. 

It  is  not  clear  at  tliis  time  whether  or  not  the  fast  method,  Algorithm  2,  is  useful  in  a  practi¬ 
cal  environment.  Initial  experiments  performed  by  Verheijen  [VS? 83]  indicate  that  the  fast 
method,  Algorithm  2,  outperforms  the  classical  one,  Algorithm  3,  only  when  n  is  greater  than 
approximately  1600.  However,  as  for  the  fast  methods  for  polynomial  multiplication,  a  hybrid  of 
Algorithm  2  and  Algorithm  3  should  significantly  lower  this  cross-over  point.  This  remains  a  sub¬ 
ject  for  future  research. 
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An  effective  implementation  of  an  off-diagonal  algorithm  should  prove  to  be  a  substantial 
tool  in  the  design  of  a  symbolic  and  algebraic  system.  A  single  routine  serves 

(1)  to  obtain  rational  approx imants  of  power  series, 

(2)  to  convert  rational  funtions  from  their  power  series  representation, 

(3)  to  compute  greatest  common  divisors  of  polynomials,  and 

(4)  to  solve  Hankel  and  Toeplitz  systems  of  equations. 

The  scope  of  this  thesis  includes  those  power  series  with  coefficients  over  a  field,  only.  This 
restriction  is  relevant  whenever  division  is  required  by  the  off-diagonal  algorithms.  By  choosing  to 
perform  pseudo-division,  rather  than  division,  the  algorithms  can  be  modified  so  that  they  are 
applicable  to  power  series  over  a  Euclidean  domain,  rather  than  a  field.  As  for  Euclid’s  extended 
algorithm,  the  modified  algorithms  shall  experience  exponential  growth  of  coefficients.  It  is  a  sub¬ 
ject  for  future  research  to  determine  if  methods  similar  to  those  of  Collins  [COLGl]  can  be  used  to 
keep  the  growth  linear.  Other  plausible  methods  for  extending  the  results  to  Euclidean  domains, 
and  which  need  to  be  investigated,  include  the  Chinese  Remainder  and  Hensel  algorithms. 

The  extension  of  results  to  Euclidean  domains  includes  as  a  special  case  the  multivariate 
power  series.  It  is  of  interest  to  determine  how  this  extension  would  compare  with  current 
methods  for  solving  block  Hankel  systems,  and  for  constructing  Pade  approx  imants  for  multi  variate 
power  series.  This  remains  a  subject  for  future  research. 

As  a  final  suggestion  for  future  research,  the  stability  of  the  algorithms  for  the  numerical 
computation  of  scaled  Pade  fractions  requires  investigation.  Some  very  preliminary  results  are 


encouraging. 
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